﻿#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "EasyXMen_sm2_InternalHandle.h"

/*------------------------------Internal Macro declaration.------------------------------*/
#define xmen_abs(x) ((x) < 0 ? (-(x)) : (x))
/* If the number of digits exceeds this number, the square operation will be faster. */
#define SQR_FASTER_THRESHOLD 5

/*------------------------------Internal function declaration.------------------------------*/
static XMen_SM2Type *get_mip();
static void convert(int n, XMen_BigTPtr x);
static flash mirvar_mem_variable(char *mem, int index, int sz);
static flash mirvar_mem(char *mem, int index);
static int EasyXMen_sm2_size(XMen_BigTPtr x);
static void nres_premult(XMen_BigTPtr x, int k, XMen_BigTPtr w);
static void nres_modadd(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr w);
static int jack(XMen_BigTPtr a, XMen_BigTPtr n);
static void denom(flash x, XMen_BigTPtr y);
static void insign(int s, flash x);
static XMen_SmallType xmen_shiftbits(XMen_SmallType x, int n);
static void fpack(XMen_BigTPtr n, XMen_BigTPtr d, flash x);
static void cbase(XMen_BigTPtr x, XMen_SmallType oldbase, XMen_BigTPtr y);
static int getdig(XMen_BigTPtr x, int i);
static void putdig(int n, XMen_BigTPtr x, int i);
static int instr(flash x, char *string);
static void xmen_lzero(XMen_BigTPtr x);
#ifdef XMEN_FP_ROUNDING
static XMen_SmallType xmen_sdiv(XMen_BigTPtr x, XMen_SmallType sn, XMen_LargeType isn, XMen_BigTPtr z);
#else
static XMen_SmallType xmen_sdiv(XMen_BigTPtr x, XMen_SmallType sn, XMen_BigTPtr z);
#endif
static XMen_SmallType normalise(XMen_BigTPtr x, XMen_BigTPtr y);
static int subdiv(XMen_BigTPtr x, int n, XMen_BigTPtr z);
static int remain(XMen_BigTPtr x, int n);
static XMen_SmallType muldvm(XMen_SmallType a, XMen_SmallType c, XMen_SmallType m, XMen_SmallType *rp);
static XMen_SmallType qdiv(XMen_LargeType u, XMen_LargeType v);
static void negify(flash x, flash y);
static int invmodp(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr z);
static void xmen_shift(XMen_BigTPtr x, int n, XMen_BigTPtr w);
static void nres(XMen_BigTPtr x, XMen_BigTPtr y);
static XMen_SmallType imuldiv(XMen_SmallType a, XMen_SmallType b, XMen_SmallType c, XMen_SmallType m, XMen_LargeType im, XMen_SmallType *rp);
static void redc(XMen_BigTPtr x, XMen_BigTPtr y);
static void nres_modmult(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr w);
static void epoint_getrhs(XMen_BigTPtr x, XMen_BigTPtr y);
static void decr(XMen_BigTPtr x, int n, XMen_BigTPtr z);
static int logb2(XMen_BigTPtr x);
static int xmen_testbit(XMen_BigTPtr x, int n);
static void nres_modsub(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr w);
static void expb2(int n, XMen_BigTPtr x);
static void nres_lucas(XMen_BigTPtr p, XMen_BigTPtr r, XMen_BigTPtr vp, XMen_BigTPtr v);
static int nres_moddiv(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr w);
static BOOL nres_sqroot(XMen_BigTPtr x, XMen_BigTPtr w);
static XMen_SmallType brand(void);
static void epoint_copy(XMen_ECurvePointT *a, XMen_ECurvePointT *b);
static XMen_ECurvePointT *epoint_init_mem_variable(char *mem, int index, int sz);
static void epoint_negate(XMen_ECurvePointT *p);
static XMen_ECurvePointT *epoint_init_mem(char *mem, int index);
static void *ecp_memalloc(int num);
static void nres_negate(XMen_BigTPtr x, XMen_BigTPtr w);
static void nres_div2(XMen_BigTPtr x, XMen_BigTPtr w);
static BOOL ecurve_padd(XMen_ECurvePointT *p, XMen_ECurvePointT *pa);
static int ecurve_add(XMen_ECurvePointT *p, XMen_ECurvePointT *pa);
static BOOL nres_multi_inverse(int m, XMen_BigTPtr *x, XMen_BigTPtr *w);
static BOOL epoint_multi_norm(int m, XMen_BigTPtr *work, XMen_ECurvePointT **p);
static int xmen_naf_window(XMen_BigTPtr x, XMen_BigTPtr x3, int i, int *nbs, int *nzs, int store);
static int ecurve_sub(XMen_ECurvePointT *p, XMen_ECurvePointT *pa);
static void ecp_memkill(char *mem, int num);
static void memkill(char *mem, int len);
static BOOL epoint_norm(XMen_ECurvePointT *p);
static void set_io_buffer_size(int len);
static void xmen_addbit(XMen_BigTPtr x, int n);
static void xmen_jsf(XMen_BigTPtr k0, XMen_BigTPtr k1, XMen_BigTPtr u0p, XMen_BigTPtr u0m, XMen_BigTPtr u1p, XMen_BigTPtr u1m);
static BOOL ecurve_add_sub(XMen_ECurvePointT *P, XMen_ECurvePointT *Q, XMen_ECurvePointT *PP, XMen_ECurvePointT *PM);
static void EasyXMen_berror(int nerr);
static XMen_SM2Type *EasyXMen_first_alloc();
static void *EasyXMen_alloc(int num, int size);
static void EasyXMen_free(void *addr);
static void *EasyXMen_memalloc(int num);
static XMen_SmallType EasyXMen_setbase(XMen_SmallType nb);
static XMen_SM2Type *EasyXMen_mirsys_basic(XMen_SM2Type *xmen_mip, int nd, XMen_SmallType nb);
static void zero(flash x);
static BOOL xmen_notint(flash x);
static void copy(flash x, flash y);
static XMen_SmallType muldvd(XMen_SmallType a, XMen_SmallType b, XMen_SmallType c, XMen_SmallType *rp);
static XMen_SmallType muldiv(XMen_SmallType a, XMen_SmallType b, XMen_SmallType c, XMen_SmallType m, XMen_SmallType *rp);
static void xmen_pmul(XMen_BigTPtr x, XMen_SmallType sn, XMen_BigTPtr z);
static void uconvert(unsigned int n, XMen_BigTPtr x);
static int exsign(flash x);
static void xmen_padd(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr z);
static void xmen_psub(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr z);
static void xmen_select(XMen_BigTPtr x, int d, XMen_BigTPtr y, XMen_BigTPtr z);
static void premult(XMen_BigTPtr x, int n, XMen_BigTPtr z);
static void numer(flash x, XMen_BigTPtr y);

/*------------------------------Internal Variable declaration.------------------------------*/
#ifdef XMEN_STATIC
static XMen_SM2Type mip;
static XMen_SM2Type *xmen_mip = &mip;
#else
static XMen_SM2Type *xmen_mip = NULL; /* only global variable */
#endif

/*------------------------------Global function define.------------------------------*/
#ifdef XMEN_GENERIC_AND_STATIC
XMen_SM2Type *XMen_Mirsys(XMen_SM2Type *xmen_mip, int nd, XMen_SmallType nb)
#else
XMen_SM2Type *XMen_Mirsys(int nd, XMen_SmallType nb)
#endif
{ /*  Initialize MIRACL system to   *
   *  use numbers to base nb, and   *
   *  nd digits or (-nd) bytes long */

    /* In these cases xmen_mip is passed as the first parameter */

#ifdef XMEN_GENERIC_AND_STATIC
    return mirsys_basic(xmen_mip, nd, nb);
#endif

    xmen_mip = EasyXMen_first_alloc();
    xmen_mip = get_mip();
    return EasyXMen_mirsys_basic(xmen_mip, nd, nb);
}

/* initialize big/flash number */
/**
 * @description Allocates and initializes a new multi-precision integer variable for cryptographic operations.
 *        This function sets up a `flash` type variable used in various big number arithmetic operations.
 *        It ensures that the system is active and error-free before proceeding with allocation.
 * @param[in] iv : The initial value to set for the newly created variable.
 * @return flash : A pointer to the newly allocated `flash` variable, or NULL if an error occurs.
 * @retval non-NULL : A valid `flash` variable was successfully allocated and initialized.
 * @retval NULL    : An error occurred during allocation or the system is not active.
 */
flash XMen_Mirvar(int iv)
{
    flash x;
    int align;
    char *ptr;

    if (xmen_mip->ERNUM)
        return NULL;
    XMEN_IN(23);

    if (!(xmen_mip->active))
    {
        EasyXMen_berror(XMEN_ERR_NO_MIRSYS);
        XMEN_OUT
        return NULL;
    }

    /* OK, now I control alignment.... */

    /* Allocate space for big, the length, the pointer, and the array */
    /* Do it all in one memory allocation - this is quicker */
    /* Ensure that the array has correct alignment */

    x = (XMen_BigTPtr)EasyXMen_alloc(xmen_size(xmen_mip->nib - 1), 1);
    if (x == NULL)
    {
        XMEN_OUT
        return x;
    }

    ptr = (char *)&x->w;
    align = (unsigned long)(ptr + sizeof(XMen_SmallType *)) % sizeof(XMen_SmallType);

    x->w = (XMen_SmallType *)(ptr + sizeof(XMen_SmallType *) + sizeof(XMen_SmallType) - align);

    if (iv != 0)
        convert(iv, x);
    XMEN_OUT
    return x;
}

/**
 * @description Returns the length of a multi-precision integer (`flash` or `big`) in words.
 *        This function calculates and returns the length of the given multi-precision integer,
 *        taking into account the specific representation used (`flash` or `big`).
 * @param[in] x : Pointer to the multi-precision integer whose length is to be determined.
 * @return int  : The length of the multi-precision integer in words.
 * @retval >0   : The valid length of the multi-precision integer.
 * @retval <=0  : Indicates an error or invalid input.
 */
int xmen_lent(flash x)
{ /* return length of big or flash in words */
    unsigned int lx;
    lx = (x->len & (XMEN_OBITS));
#ifdef XMEN_FLASH
    return (int)((lx & (XMEN_MSK)) + ((lx >> (XMEN_BTS)) & (XMEN_MSK)));
#else
    return (int)lx;
#endif
}

/**
 * @description Compares two multi-precision integers (`XMen_BigTPtr` type) and returns the result of the comparison.
 *        This function determines whether one multi-precision integer is greater than, less than, or equal to another.
 * @param[in] x : Pointer to the first multi-precision integer to compare.
 * @param[in] y : Pointer to the second multi-precision integer to compare.
 * @return int  : The result of the comparison.
 * @retval 1    : If `x` is greater than `y`.
 * @retval -1   : If `x` is less than `y`.
 * @retval 0    : If `x` and `y` are equal.
 */
int XMen_Compare(XMen_BigTPtr x, XMen_BigTPtr y)
{ /* compare x and y: =1 if x>y  =-1 if x<y *
   *  =0 if x=y                             */
    int m, n, sig;
    unsigned int sx, sy;
    if (x == y)
        return 0;
    sx = (x->len & XMEN_MSBIT);
    sy = (y->len & XMEN_MSBIT);
    if (sx == 0)
        sig = PLUS;
    else
        sig = MINUS;
    if (sx != sy)
        return sig;
    m = (int)(x->len & XMEN_OBITS);
    n = (int)(y->len & XMEN_OBITS);
    if (m > n)
        return sig;
    if (m < n)
        return -sig;
    while (m > 0)
    { /* check digit by digit */
        m--;
        if (x->w[m] > y->w[m])
            return sig;
        if (x->w[m] < y->w[m])
            return -sig;
    }
    return 0;
}

/**
 * @description Adds an integer to a multi-precision integer (`XMen_BigTPtr`) and stores the result in another multi-precision integer.
 *        This function performs the operation `z = x + n`, where `x` is a multi-precision integer, `n` is an integer,
 *        and `z` is the resulting multi-precision integer.
 * @param[in] x : Pointer to the multi-precision integer to which the integer `n` will be added.
 * @param[in] n : The integer value to add to the multi-precision integer `x`.
 * @param[out] z: Pointer to the multi-precision integer where the result of the addition will be stored.
 * @return void
 */
void XMen_Incr(XMen_BigTPtr x, int n, XMen_BigTPtr z)
{ /* XMen_Add int to big number: z=x+n */

    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(7)

    convert(n, xmen_mip->w0);
    xmen_select(x, PLUS, xmen_mip->w0, z);

    XMEN_OUT
}

/**
 * @description Converts a sequence of bytes into a multi-precision integer (`XMen_BigTPtr`).
 *        This function interprets the byte array as a big-endian representation of an integer,
 *        with the first byte being the most significant.
 * @param[in] len : The number of bytes to convert.
 * @param[in] ptr : Pointer to the byte array containing the data to convert.
 * @param[out] x  : Pointer to the multi-precision integer where the result will be stored.
 * @return void
 */
void XMen_Bytes2Big(int len, const char *ptr, XMen_BigTPtr x)
{ /* convert len bytes into a big           *
   * The first byte is the Most significant */
    int i, j, m, n, r;
    unsigned int dig;
    unsigned char ch;
    XMen_SmallType wrd;

    if (xmen_mip->ERNUM)
        return;
    XMEN_IN(140);

    zero(x);

    if (len <= 0)
    {
        XMEN_OUT
        return;
    }
    /* remove leading zeros.. */

    while (*ptr == 0)
    {
        ptr++;
        len--;
        if (len == 0)
        {
            XMEN_OUT
            return;
        }
    }

#ifndef XMEN_SIMPLE_BASE
    if (xmen_mip->base == 0)
    { /* pack bytes directly into big */
#endif
#ifndef XMEN_NOFULLWIDTH
        m = MIRACL / 8;
        n = len / m;

        r = len % m;
        wrd = (XMen_SmallType)0;
        if (r != 0)
        {
            n++;
            for (j = 0; j < r; j++)
            {
                wrd <<= 8;
                wrd |= XMEN_TOBYTE(*ptr++);
            }
        }
        x->len = n;
        if (n > xmen_mip->nib && xmen_mip->check)
        {
            EasyXMen_berror(XMEN_ERR_OVERFLOW);
            XMEN_OUT
            return;
        }
        if (r != 0)
        {
            n--;
            x->w[n] = wrd;
        }

        for (i = n - 1; i >= 0; i--)
        {
            for (j = 0; j < m; j++)
            {
                wrd <<= 8;
                wrd |= XMEN_TOBYTE(*ptr++);
            }
            x->w[i] = wrd;
        }
        xmen_lzero(x); /* needed */
#endif
#ifndef XMEN_SIMPLE_BASE
    }
    else
    {
        for (i = 0; i < len; i++)
        {
            if (xmen_mip->ERNUM)
                break;
#if MIRACL == 8
            xmen_shift(x, 1, x);
#else
            premult(x, 256, x);
#endif
            ch = XMEN_TOBYTE(ptr[i]);
            dig = ch;
            XMen_Incr(x, (int)dig, x);
        }
    }
#endif
    XMEN_OUT
}

/**
 * @description Adds two signed multi-precision integers (`XMen_BigTPtr`) together and stores the result in a third multi-precision integer.
 *        This function performs the operation `z = x + y`, where `x` and `y` are multi-precision integers,
 *        and `z` is the resulting multi-precision integer.
 * @param[in] x : Pointer to the first multi-precision integer to add.
 * @param[in] y : Pointer to the second multi-precision integer to add.
 * @param[out] z: Pointer to the multi-precision integer where the result of the addition will be stored.
 * @return void
 */
void XMen_Add(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr z)
{ /* XMen_Add two signed big numbers together z=x+y */

    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(27)
    xmen_select(x, PLUS, y, z);
    XMEN_OUT
}

/**
 * @description Converts a string representation of a number in a specified base (`IOBASE`) into a multi-precision integer (`flash`).
 *        This function parses the input string and stores the resulting multi-precision integer in the provided variable.
 * @param[in,out] x     : Pointer to the multi-precision integer where the result will be stored.
 * @param[in]     string: Null-terminated string containing the number to convert, represented in base `IOBASE`.
 * @return int         : Returns a status code indicating success or failure.
 * @retval 0           : An error occurred during conversion or the system is not ready.
 * @retval non-zero    : Conversion was successful.
 */
int XMen_Cinstr(flash x, char *string)
{ /*  input big number in base IOBASE  */
    XMen_SmallType newb, oldb, b;
    unsigned int lx;
    int ipt;

    if (xmen_mip->ERNUM)
        return 0;

    XMEN_IN(78)

    newb = xmen_mip->IOBASE;
    oldb = xmen_mip->apbase;
    EasyXMen_setbase(newb); /* temporarily change base ... */
    b = xmen_mip->base;
    xmen_mip->check = OFF;
    ipt = instr(xmen_mip->w5, string); /* ... and get number  */

    xmen_mip->check = ON;
    lx = (xmen_mip->w5->len & XMEN_OBITS);
#ifdef XMEN_FLASH
    if ((int)(lx & XMEN_MSK) > xmen_mip->nib || (int)((lx >> XMEN_BTS) & XMEN_MSK) > xmen_mip->nib)
#else
    if ((int)lx > xmen_mip->nib)
#endif
    { /* numerator or denominator too big */
        EasyXMen_berror(XMEN_ERR_OVERFLOW);
        XMEN_OUT
        return 0;
    }
    EasyXMen_setbase(oldb); /* restore original base */

    cbase(xmen_mip->w5, b, x);

    XMEN_OUT
    return ipt;
}

/**
 * @description Divides two multi-precision integers (`XMen_BigTPtr`) and optionally returns the quotient or remainder.
 *        This function performs the operation `z = x / y` for the quotient, and `x = x % y` for the remainder,
 *        depending on the parameters passed. It can return either the quotient or the remainder based on the
 *        provided pointers.
 * @param[in,out] x : Pointer to the first multi-precision integer (dividend). After the operation, if `y` is not `z`,
 *                    `x` will contain the remainder of the division (`x = x % y`).
 * @param[in]     y : Pointer to the second multi-precision integer (divisor).
 * @param[out]    z : Pointer to the multi-precision integer where the result (quotient) will be stored. If `z` is
 *                    equal to `x`, it will store the quotient; if `z` is equal to `y`, it will store the remainder.
 * @return void
 */
void XMen_Divide(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr z)
{ /*  divide two big numbers  z=x/y : x=x mod y  *
   *  returns quotient only if  XMen_Divide(x,y,x)    *
   *  returns remainder only if XMen_Divide(x,y,y)    */
    XMen_SmallType carry, attemp, ldy, sdy, ra, r, d, tst, psum;
    unsigned int sx, sy, sz;
    XMen_SmallType borrow, dig, *w0g, *yg;
    int i, k, m, x0, y0, w00;
    XMen_BigTPtr w0;

#ifdef XMEN_ITANIUM
    XMen_SmallType tm;
#endif
#ifdef XMEN_WIN64
    XMen_SmallType tm;
#endif
    union doubleword dble;
    XMen_LargeType dbled;
    XMen_LargeType ldres;
    BOOL check;

    if (xmen_mip->ERNUM)
        return;
    w0 = xmen_mip->w0;

    XMEN_IN(6)

    if (x == y)
        EasyXMen_berror(XMEN_ERR_BAD_PARAMETERS);
#ifdef XMEN_FLASH
    if (xmen_notint(x) || xmen_notint(y))
        EasyXMen_berror(XMEN_ERR_INT_OP);
#endif
    if (y->len == 0)
        EasyXMen_berror(XMEN_ERR_DIV_BY_ZERO);
    if (xmen_mip->ERNUM)
    {
        XMEN_OUT
        return;
    }
    sx = (x->len & XMEN_MSBIT); /*  extract signs ... */
    sy = (y->len & XMEN_MSBIT);
    sz = (sx ^ sy);
    x->len &= XMEN_OBITS; /* ... and force operands to positive  */
    y->len &= XMEN_OBITS;
    x0 = (int)x->len;
    y0 = (int)y->len;
    copy(x, w0);
    w00 = (int)w0->len;
    if (xmen_mip->check && (w00 - y0 + 1 > xmen_mip->nib))
    {
        EasyXMen_berror(XMEN_ERR_OVERFLOW);
        XMEN_OUT
        return;
    }
    d = 0;
    if (x0 == y0)
    {
        if (x0 == 1) /* special case - x and y are both xmen_smalls */
        {
            d = XMEN_DIV(w0->w[0], y->w[0]);
            w0->w[0] = XMEN_REMAIN(w0->w[0], y->w[0]);
            xmen_lzero(w0);
        }
        else if (XMEN_DIV(w0->w[x0 - 1], 4) < y->w[x0 - 1])
            while (XMen_Compare(w0, y) >= 0)
            { /* XMen_SmallType quotient - so do up to four subtracts instead */
                xmen_psub(w0, y, w0);
                d++;
            }
    }
    if (XMen_Compare(w0, y) < 0)
    {               /*  x less than y - so x becomes remainder */
        if (x != z) /* testing parameters */
        {
            copy(w0, x);
            if (x->len != 0)
                x->len |= sx;
        }
        if (y != z)
        {
            zero(z);
            z->w[0] = d;
            if (d > 0)
                z->len = (sz | 1);
        }
        y->len |= sy;
        XMEN_OUT
        return;
    }

    if (y0 == 1)
    { /* y is int - so use subdiv instead */
#ifdef XMEN_FP_ROUNDING
        r = xmen_sdiv(w0, y->w[0], xmen_invert(y->w[0]), w0);
#else
        r = xmen_sdiv(w0, y->w[0], w0);
#endif
        if (y != z)
        {
            copy(w0, z);
            z->len |= sz;
        }
        if (x != z)
        {
            zero(x);
            x->w[0] = r;
            if (r > 0)
                x->len = (sx | 1);
        }
        y->len |= sy;
        XMEN_OUT
        return;
    }
    if (y != z)
        zero(z);
    d = normalise(y, y);
    check = xmen_mip->check;
    xmen_mip->check = OFF;
#ifndef XMEN_SIMPLE_BASE
    if (xmen_mip->base == 0)
    {
#endif
#ifndef XMEN_NOFULLWIDTH
        if (d != 1)
            xmen_pmul(w0, d, w0);
        ldy = y->w[y0 - 1];
        sdy = y->w[y0 - 2];
        w0g = w0->w;
        yg = y->w;
        for (k = w00 - 1; k >= y0 - 1; k--)
        { /* long division */

            carry = 0;
            if (w0->w[k + 1] == ldy) /* guess next quotient digit */
            {
                attemp = (XMen_SmallType)(-1);
                ra = ldy + w0->w[k];
                if (ra < ldy)
                    carry = 1;
            }
            else
            {
                dble.h[BYTE_BOT] = w0->w[k];
                dble.h[BYTE_TOP] = w0->w[k + 1];
                attemp = (XMen_SmallType)(dble.d / ldy);
                ra = (XMen_SmallType)(dble.d - (XMen_LargeType)attemp * ldy);
            }
            while (carry == 0)
            {
                dble.d = (XMen_LargeType)attemp * sdy;
                r = dble.h[BYTE_BOT];
                tst = dble.h[BYTE_TOP];
                if (tst < ra || (tst == ra && r <= w0->w[k - 1]))
                    break;
                attemp--; /* refine guess */
                ra += ldy;
                if (ra < ldy)
                    carry = 1;
            }
            m = k - y0 + 1;
            if (attemp > 0)
            { /* do partial subtraction */
                borrow = 0;
                /*  inline - substitutes for loop below */

                for (i = 0; i < y0; i++)
                {
                    dble.d = (XMen_LargeType)attemp * y->w[i] + borrow;
                    dig = dble.h[BYTE_BOT];
                    borrow = dble.h[BYTE_TOP];
                    if (w0->w[m + i] < dig)
                        borrow++;
                    w0->w[m + i] -= dig;
                }

                if (w0->w[k + 1] < borrow)
                { /* whoops! - over did it */
                    w0->w[k + 1] = 0;
                    carry = 0;
                    for (i = 0; i < y0; i++)
                    { /* compensate for error ... */
                        psum = w0->w[m + i] + y->w[i] + carry;
                        if (psum > y->w[i])
                            carry = 0;
                        if (psum < y->w[i])
                            carry = 1;
                        w0->w[m + i] = psum;
                    }
                    attemp--; /* ... and adjust guess */
                }
                else
                    w0->w[k + 1] -= borrow;
            }
            if (k == w00 - 1 && attemp == 0)
                w00--;
            else if (y != z)
                z->w[m] = attemp;
        }
#endif
#ifndef XMEN_SIMPLE_BASE
    }
    else
    { /* have to do it the hard way */
        if (d != 1)
            xmen_pmul(w0, d, w0);
        ldy = y->w[y0 - 1];
        sdy = y->w[y0 - 2];

        for (k = w00 - 1; k >= y0 - 1; k--)
        { /* long division */

            if (w0->w[k + 1] == ldy) /* guess next quotient digit */
            {
                attemp = xmen_mip->base - 1;
                ra = ldy + w0->w[k];
            }
            else
            {
                dbled = (XMen_LargeType)w0->w[k + 1] * xmen_mip->base + w0->w[k];
                attemp = (XMen_SmallType)XMEN_LROUND(dbled / ldy);
                ra = (XMen_SmallType)(dbled - (XMen_LargeType)attemp * ldy);
            }
            while (ra < xmen_mip->base)
            {
                dbled = (XMen_LargeType)sdy * attemp;
#ifdef XMEN_FP_ROUNDING
                tst = (XMen_SmallType)XMEN_LROUND(dbled * xmen_mip->inverse_base);
#else
                if (xmen_mip->base == xmen_mip->base2)
                    tst = (XMen_SmallType)(dbled >> xmen_mip->lg2b);
                else
                    tst = (XMen_SmallType)XMEN_LROUND(dbled / xmen_mip->base);
#endif
                r = (XMen_SmallType)(dbled - (XMen_LargeType)tst * xmen_mip->base);

                if (tst < ra || (tst == ra && r <= w0->w[k - 1]))
                    break;
                attemp--; /* refine guess */
                ra += ldy;
            }
            m = k - y0 + 1;
            if (attemp > 0)
            { /* do partial subtraction */
                borrow = 0;
                for (i = 0; i < y0; i++)
                {
                    dbled = (XMen_LargeType)attemp * y->w[i] + borrow;
#ifdef XMEN_FP_ROUNDING
                    borrow = (XMen_SmallType)XMEN_LROUND(dbled * xmen_mip->inverse_base);
#else
                    if (xmen_mip->base == xmen_mip->base2)
                        borrow = (XMen_SmallType)(dbled >> xmen_mip->lg2b);
                    else
                        borrow = (XMen_SmallType)XMEN_LROUND(dbled / xmen_mip->base);
#endif
                    dig = (XMen_SmallType)(dbled - (XMen_LargeType)borrow * xmen_mip->base);

                    if (w0->w[m + i] < dig)
                    { /* set borrow */
                        borrow++;
                        w0->w[m + i] += (xmen_mip->base - dig);
                    }
                    else
                        w0->w[m + i] -= dig;
                }
                if (w0->w[k + 1] < borrow)
                { /* whoops! - over did it */
                    w0->w[k + 1] = 0;
                    carry = 0;
                    for (i = 0; i < y0; i++)
                    { /* compensate for error ... */
                        psum = w0->w[m + i] + y->w[i] + carry;
                        carry = 0;
                        if (psum >= xmen_mip->base)
                        {
                            carry = 1;
                            psum -= xmen_mip->base;
                        }
                        w0->w[m + i] = psum;
                    }
                    attemp--; /* ... and adjust guess */
                }
                else
                    w0->w[k + 1] -= borrow;
            }
            if (k == w00 - 1 && attemp == 0)
                w00--;
            else if (y != z)
                z->w[m] = attemp;
        }
    }
#endif
    if (y != z)
        z->len = ((w00 - y0 + 1) | sz); /* set sign and length of result */

    w0->len = y0;

    xmen_lzero(y);
    xmen_lzero(z);

    if (x != z)
    {
        xmen_lzero(w0);
#ifdef XMEN_FP_ROUNDING
        if (d != 1)
            xmen_sdiv(w0, d, xmen_invert(d), x);
#else
        if (d != 1)
            xmen_sdiv(w0, d, x);
#endif
        else
            copy(w0, x);
        if (x->len != 0)
            x->len |= sx;
    }
#ifdef XMEN_FP_ROUNDING
    if (d != 1)
        xmen_sdiv(y, d, xmen_invert(d), y);
#else
    if (d != 1)
        xmen_sdiv(y, d, y);
#endif
    y->len |= sy;
    xmen_mip->check = check;

    XMEN_OUT
}

/**
 * @description Subtracts one multi-precision integer from another and stores the result in a third multi-precision integer.
 *        This function performs the operation `z = x - y`, where `x` and `y` are multi-precision integers,
 *        and `z` is the resulting multi-precision integer.
 * @param[in] x : Pointer to the first multi-precision integer (minuend) from which the second number will be subtracted.
 * @param[in] y : Pointer to the second multi-precision integer (subtrahend) to subtract from the first number.
 * @param[out] z: Pointer to the multi-precision integer where the result of the subtraction will be stored.
 * @return void
 */
void XMen_Subtract(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr z)
{ /* XMen_Subtract two big signed numbers z=x-y */

    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(28)

    xmen_select(x, MINUS, y, z);

    XMEN_OUT
}

/**
 * @description Performs a combined operation of multiplication, addition, and division on multi-precision integers.
 *        This function calculates `q = (x * y + z) / w` and optionally returns the remainder in `r`.
 *        It can return either the quotient or the remainder based on the provided pointers.
 *        Addition is performed only if `x`, `y`, and `z` are distinct.
 * @param[in] x : Pointer to the first multi-precision integer to multiply.
 * @param[in] y : Pointer to the second multi-precision integer to multiply.
 * @param[in] z : Pointer to the multi-precision integer to add after multiplication.
 * @param[in] w : Pointer to the multi-precision integer divisor for the final division.
 * @param[out] q: Pointer to the multi-precision integer where the quotient will be stored.
 *                If `q` equals `r`, only the remainder is returned.
 * @param[out] r: Pointer to the multi-precision integer where the remainder will be stored.
 *                If `q` equals `r`, only the quotient is returned.
 * @return void
 */
void XMen_Mad(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr z, XMen_BigTPtr w, XMen_BigTPtr q, XMen_BigTPtr r)
{ /* Multiply, Add and Divide; q=(x*y+z)/w remainder r   *
   * returns remainder only if w=q, quotient only if q=r *
   * XMen_Add done only if x, y and z are distinct.           */

    BOOL check;
    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(24)
    if (w == r)
    {
        EasyXMen_berror(XMEN_ERR_BAD_PARAMETERS);
        XMEN_OUT
        return;
    }
    check = xmen_mip->check;
    xmen_mip->check = OFF; /* turn off some error checks */

    XMen_Multiply(x, y, xmen_mip->w0);
    if (x != z && y != z)
        XMen_Add(xmen_mip->w0, z, xmen_mip->w0);

    XMen_Divide(xmen_mip->w0, w, q);
    if (q != r)
        copy(xmen_mip->w0, r);
    xmen_mip->check = check;
    XMEN_OUT
}

/**
 * @description Calculates the greatest common divisor (GCD) of two multi-precision integers using Euclid's method.
 *        Additionally, it calculates coefficients `xd` and `yd` such that `z = x*xd + y*yd = gcd(x,y)`.
 *        The function returns the GCD in `z` if `z` is distinct from `xd` and `yd`, and ensures that `xd` is positive
 *        while `yd` is negative. If `xd` and `yd` are not distinct, only `xd` is calculated.
 * @param[in] x   : Pointer to the first multi-precision integer.
 * @param[in] y   : Pointer to the second multi-precision integer.
 * @param[out] xd : Pointer to the multi-precision integer where the coefficient for `x` will be stored.
 * @param[out] yd : Pointer to the multi-precision integer where the coefficient for `y` will be stored.
 *                  If `xd` and `yd` are the same, only `xd` is calculated.
 * @param[out] z  : Pointer to the multi-precision integer where the GCD will be stored.
 *                  Only returned if distinct from `xd` and `yd`.
 * @return int    : Returns the size of the result (`z`) or an error code.
 * @retval 0      : An error occurred during the calculation or system is not ready.
 * @retval non-zero: Indicates the size of the GCD result.
 */
int XMen_Xgcd(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr xd, XMen_BigTPtr yd, XMen_BigTPtr z)
{ /* greatest common divisor by Euclids method  *
   * extended to also calculate xd and yd where *
   *      z = x.xd + y.yd = gcd(x,y)            *
   * if xd, yd not distinct, only xd calculated *
   * z only returned if distinct from xd and yd *
   * xd will always be positive, yd negative    */

    int s, n, iter;
    XMen_SmallType r, a, b, c, d;
    XMen_SmallType q, m, sr;

#ifdef XMen_DltType
    union doubleword uu, vv;
    XMen_LargeType u, v, lr;
#else
    XMen_SmallType u, v, lr;
#endif

    BOOL last, dplus = TRUE;
    XMen_BigTPtr t;

    if (xmen_mip->ERNUM)
        return 0;

    XMEN_IN(30)

#ifdef XMEN_COUNT_OPS
    fpx++;
#endif

    copy(x, xmen_mip->w1);
    copy(y, xmen_mip->w2);
    s = exsign(xmen_mip->w1);
    insign(PLUS, xmen_mip->w1);
    insign(PLUS, xmen_mip->w2);
    convert(1, xmen_mip->w3);
    zero(xmen_mip->w4);
    last = FALSE;
    a = b = c = d = 0;
    iter = 0;

    while (EasyXMen_sm2_size(xmen_mip->w2) != 0)
    {
        if (b == 0)
        { /* update xmen_mip->w1 and xmen_mip->w2 */

            XMen_Divide(xmen_mip->w1, xmen_mip->w2, xmen_mip->w5);
            t = xmen_mip->w1, xmen_mip->w1 = xmen_mip->w2, xmen_mip->w2 = t; /* swap(xmen_mip->w1,xmen_mip->w2) */
            XMen_Multiply(xmen_mip->w4, xmen_mip->w5, xmen_mip->w0);
            XMen_Add(xmen_mip->w3, xmen_mip->w0, xmen_mip->w3);
            t = xmen_mip->w3, xmen_mip->w3 = xmen_mip->w4, xmen_mip->w4 = t; /* swap(xd,yd) */
            iter++;
        }
        else
        {

            /* printf("a= %I64u b= %I64u c= %I64u  d= %I64u \n",a,b,c,d);   */

            xmen_pmul(xmen_mip->w1, c, xmen_mip->w5); /* c*w1 */
            xmen_pmul(xmen_mip->w1, a, xmen_mip->w1); /* a*w1 */
            xmen_pmul(xmen_mip->w2, b, xmen_mip->w0); /* b*w2 */
            xmen_pmul(xmen_mip->w2, d, xmen_mip->w2); /* d*w2 */

            if (!dplus)
            {
                xmen_psub(xmen_mip->w0, xmen_mip->w1, xmen_mip->w1); /* b*w2-a*w1 */
                xmen_psub(xmen_mip->w5, xmen_mip->w2, xmen_mip->w2); /* c*w1-d*w2 */
            }
            else
            {
                xmen_psub(xmen_mip->w1, xmen_mip->w0, xmen_mip->w1); /* a*w1-b*w2 */
                xmen_psub(xmen_mip->w2, xmen_mip->w5, xmen_mip->w2); /* d*w2-c*w1 */
            }
            xmen_pmul(xmen_mip->w3, c, xmen_mip->w5);
            xmen_pmul(xmen_mip->w3, a, xmen_mip->w3);
            xmen_pmul(xmen_mip->w4, b, xmen_mip->w0);
            xmen_pmul(xmen_mip->w4, d, xmen_mip->w4);

            if (a == 0)
                copy(xmen_mip->w0, xmen_mip->w3);
            else
                xmen_padd(xmen_mip->w3, xmen_mip->w0, xmen_mip->w3);
            xmen_padd(xmen_mip->w4, xmen_mip->w5, xmen_mip->w4);
        }
        if (xmen_mip->ERNUM || EasyXMen_sm2_size(xmen_mip->w2) == 0)
            break;

        n = (int)xmen_mip->w1->len;
        if (n == 1)
        {
            last = TRUE;
            u = xmen_mip->w1->w[0];
            v = xmen_mip->w2->w[0];
        }
        else
        {
            m = xmen_mip->w1->w[n - 1] + 1;
#ifndef XMEN_SIMPLE_BASE
            if (xmen_mip->base == 0)
            {
#endif
#ifndef XMEN_NOFULLWIDTH
#ifdef XMen_DltType
                /* use double length type if available */
                if (n > 2 && m != 0)
                { /* squeeze out as much significance as possible */
                    uu.h[BYTE_TOP] = muldvm(xmen_mip->w1->w[n - 1], xmen_mip->w1->w[n - 2], m, &sr);
                    uu.h[BYTE_BOT] = muldvm(sr, xmen_mip->w1->w[n - 3], m, &sr);
                    vv.h[BYTE_TOP] = muldvm(xmen_mip->w2->w[n - 1], xmen_mip->w2->w[n - 2], m, &sr);
                    vv.h[BYTE_BOT] = muldvm(sr, xmen_mip->w2->w[n - 3], m, &sr);
                }
                else
                {
                    uu.h[BYTE_TOP] = xmen_mip->w1->w[n - 1];
                    uu.h[BYTE_BOT] = xmen_mip->w1->w[n - 2];
                    vv.h[BYTE_TOP] = xmen_mip->w2->w[n - 1];
                    vv.h[BYTE_BOT] = xmen_mip->w2->w[n - 2];
                    if (n == 2)
                        last = TRUE;
                }

                u = uu.d;
                v = vv.d;
#else
                if (m == 0)
                {
                    u = xmen_mip->w1->w[n - 1];
                    v = xmen_mip->w2->w[n - 1];
                }
                else
                {
                    u = muldvm(xmen_mip->w1->w[n - 1], xmen_mip->w1->w[n - 2], m, &sr);
                    v = muldvm(xmen_mip->w2->w[n - 1], xmen_mip->w2->w[n - 2], m, &sr);
                }
#endif
#endif
#ifndef XMEN_SIMPLE_BASE
            }
            else
            {
#ifdef XMen_DltType
                if (n > 2)
                { /* squeeze out as much significance as possible */
                    u = muldiv(xmen_mip->w1->w[n - 1], xmen_mip->base, xmen_mip->w1->w[n - 2], m, &sr);
                    u = u * xmen_mip->base + muldiv(sr, xmen_mip->base, xmen_mip->w1->w[n - 3], m, &sr);
                    v = muldiv(xmen_mip->w2->w[n - 1], xmen_mip->base, xmen_mip->w2->w[n - 2], m, &sr);
                    v = v * xmen_mip->base + muldiv(sr, xmen_mip->base, xmen_mip->w2->w[n - 3], m, &sr);
                }
                else
                {
                    u = (XMen_LargeType)xmen_mip->base * xmen_mip->w1->w[n - 1] + xmen_mip->w1->w[n - 2];
                    v = (XMen_LargeType)xmen_mip->base * xmen_mip->w2->w[n - 1] + xmen_mip->w2->w[n - 2];
                    last = TRUE;
                }
#else
                u = muldiv(xmen_mip->w1->w[n - 1], xmen_mip->base, xmen_mip->w1->w[n - 2], m, &sr);
                v = muldiv(xmen_mip->w2->w[n - 1], xmen_mip->base, xmen_mip->w2->w[n - 2], m, &sr);
#endif
            }
#endif
        }

        dplus = TRUE;
        a = 1;
        b = 0;
        c = 0;
        d = 1;

        for (;;)
        { /* work only with most significant piece */
            if (last)
            {
                if (v == 0)
                    break;
                q = qdiv(u, v);
                if (q == 0)
                    break;
            }
            else
            {
                if (dplus)
                {
                    if ((XMen_SmallType)(v - c) == 0 || (XMen_SmallType)(v + d) == 0)
                        break;

                    q = qdiv(u + a, v - c);

                    if (q == 0)
                        break;

                    if (q != qdiv(u - b, v + d))
                        break;
                }
                else
                {
                    if ((XMen_SmallType)(v + c) == 0 || (XMen_SmallType)(v - d) == 0)
                        break;
                    q = qdiv(u - a, v + c);
                    if (q == 0)
                        break;
                    if (q != qdiv(u + b, v - d))
                        break;
                }
            }

            if (q == 1)
            {
                if ((XMen_SmallType)(b + d) >= MAXBASE)
                    break;
                r = a + c;
                a = c;
                c = r;
                r = b + d;
                b = d;
                d = r;
                lr = u - v;
                u = v;
                v = lr;
            }
            else
            {
                if (q >= XMEN_DIV(MAXBASE - b, d))
                    break;
                r = a + q * c;
                a = c;
                c = r;
                r = b + q * d;
                b = d;
                d = r;
                lr = u - q * v;
                u = v;
                v = lr;
            }
            iter++;
            dplus = !dplus;
        }
        iter %= 2;
    }

    if (s == MINUS)
        iter++;
    if (iter % 2 == 1)
        XMen_Subtract(y, xmen_mip->w3, xmen_mip->w3);

    if (xd != yd)
    {
        negify(x, xmen_mip->w2);
        XMen_Mad(xmen_mip->w2, xmen_mip->w3, xmen_mip->w1, y, xmen_mip->w4, xmen_mip->w4);
        copy(xmen_mip->w4, yd);
    }
    copy(xmen_mip->w3, xd);
    if (z != xd && z != yd)
        copy(xmen_mip->w1, z);

    XMEN_OUT
    return (EasyXMen_sm2_size(xmen_mip->w1));
}

/**
 * @description Initializes an elliptic curve for cryptographic operations.
 *        This function sets up the parameters `a`, `b`, and `p` of the elliptic curve equation `y^2 = x^3 + ax + b (mod p)`.
 *        It also prepares the Montgomery reduction parameters for efficient modular arithmetic and sets the coordinate system
 *        to be used for point representation on the curve. The sizes of `A` and `B` are determined and adjusted if they are
 *        too large by checking if they represent negative small numbers in the field defined by `p`.
 * @param[in] a  : Pointer to the multi-precision integer representing the coefficient `a` in the elliptic curve equation.
 * @param[in] b  : Pointer to the multi-precision integer representing the coefficient `b` in the elliptic curve equation.
 * @param[in] p  : Pointer to the multi-precision integer representing the prime modulus `p` of the finite field.
 * @param[in] type: An integer specifying the type of coordinates to use for points on the curve. If set to `XMEN_BEST`,
 *                  the function will choose the most efficient coordinate system available, otherwise it will use the specified type.
 * @return void
 */
void XMen_EcurveInit(XMen_BigTPtr a, XMen_BigTPtr b, XMen_BigTPtr p, int type)
{ /* Initialize the active ecurve    *
   * Asize indicate size of A        *
   * Bsize indicate size of B        */
    int as;

    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(93)

#ifndef XMEN_NO_SS
    xmen_mip->SS = FALSE; /* no special support for super-singular curves */
#endif

    prepare_monty(p);

    xmen_mip->Asize = EasyXMen_sm2_size(a);
    if (xmen_abs(xmen_mip->Asize) == XMEN_TOOBIG)
    {
        if (xmen_mip->Asize >= 0)
        { /* big positive number - check it isn't minus something small */
            copy(a, xmen_mip->w1);
            XMen_Divide(xmen_mip->w1, p, p);
            XMen_Subtract(p, xmen_mip->w1, xmen_mip->w1);
            as = EasyXMen_sm2_size(xmen_mip->w1);
            if (as < XMEN_TOOBIG)
                xmen_mip->Asize = -as;
        }
    }
    nres(a, xmen_mip->A);

    xmen_mip->Bsize = EasyXMen_sm2_size(b);
    if (xmen_abs(xmen_mip->Bsize) == XMEN_TOOBIG)
    {
        if (xmen_mip->Bsize >= 0)
        { /* big positive number - check it isn't minus something small */
            copy(b, xmen_mip->w1);
            XMen_Divide(xmen_mip->w1, p, p);
            XMen_Subtract(p, xmen_mip->w1, xmen_mip->w1);
            as = EasyXMen_sm2_size(xmen_mip->w1);
            if (as < XMEN_TOOBIG)
                xmen_mip->Bsize = -as;
        }
    }

    nres(b, xmen_mip->B);
#ifdef XMEN_EDWARDS
    xmen_mip->coord = XMEN_PROJECTIVE; /* only type supported for Edwards curves */
#else
    if (type == XMEN_BEST)
        xmen_mip->coord = XMEN_PROJECTIVE;
    else
        xmen_mip->coord = type;
#endif
    XMEN_OUT
    return;
}

/**
 * @description Initializes an elliptic curve point structure to represent the point at infinity.
 *        This function allocates memory for a new `XMen_ECurvePointT` structure and initializes it
 *        to represent the general point at infinity on an elliptic curve. The point at infinity is
 *        a special point that serves as the identity element for the group operation on the curve.
 * @return XMen_ECurvePointT* : Pointer to the initialized elliptic curve point structure, or NULL if an error occurs.
 * @retval NULL : An error occurred during initialization or system is not ready.
 * @retval non-NULL: A pointer to the newly allocated and initialized elliptic curve point structure.
 */
XMen_ECurvePointT *XMen_EpointInit(void)
{ /* initialise XMen_ECurvePointT to general point at infinity. */
    XMen_ECurvePointT *p;
    char *ptr;

    if (xmen_mip->ERNUM)
        return NULL;

    XMEN_IN(96)

    /* Create space for whole structure in one heap access */

    p = (XMen_ECurvePointT *)EasyXMen_alloc(xmen_esize(xmen_mip->nib - 1), 1);

    ptr = (char *)p + sizeof(XMen_ECurvePointT);
    p->X = mirvar_mem(ptr, 0);
    p->Y = mirvar_mem(ptr, 1);
    p->Z = mirvar_mem(ptr, 2);
    p->marker = XMEN_EPOINT_INFINITY;

    XMEN_OUT

    return p;
}

/**
 * @description Multiplies two multi-precision integers and stores the result in a third multi-precision integer.
 *        This function performs the operation `z = x * y`, where `x` and `y` are multi-precision integers,
 *        and `z` is the resulting multi-precision integer. It handles large number multiplication efficiently
 *        using long multiplication or optimized squaring for equal operands, ensuring that overflow does not occur.
 * @param[in] x : Pointer to the first multi-precision integer to multiply.
 * @param[in] y : Pointer to the second multi-precision integer to multiply.
 * @param[out] z: Pointer to the multi-precision integer where the result of the multiplication will be stored.
 * @return void
 */
void XMen_Multiply(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr z)
{ /*  XMen_Multiply two big numbers: z=x.y  */
    int i, xl, yl, j, ti;
    XMen_SmallType carry, *xg, *yg, *w0g;

#ifdef XMEN_ITANIUM
    XMen_SmallType tm;
#endif
#ifdef XMEN_WIN64
    XMen_SmallType tm, tr;
#endif
    unsigned int sz;
    XMen_BigTPtr w0;
    union doubleword dble;
    XMen_LargeType dbled;
    XMen_LargeType ldres;

    if (xmen_mip->ERNUM)
        return;
    if (y->len == 0 || x->len == 0)
    {
        zero(z);
        return;
    }
    if (x != xmen_mip->w5 && y != xmen_mip->w5 && z == xmen_mip->w5)
        w0 = xmen_mip->w5;
    else
        w0 = xmen_mip->w0; /* local pointer */

    XMEN_IN(5)

#ifdef XMEN_FLASH
    if (xmen_notint(x) || xmen_notint(y))
    {
        EasyXMen_berror(XMEN_ERR_INT_OP);
        XMEN_OUT
        return;
    }
#endif
    sz = ((x->len & XMEN_MSBIT) ^ (y->len & XMEN_MSBIT));
    xl = (int)(x->len & XMEN_OBITS);
    yl = (int)(y->len & XMEN_OBITS);
    zero(w0);
    if (xmen_mip->check && xl + yl > xmen_mip->nib)
    {
        EasyXMen_berror(XMEN_ERR_OVERFLOW);
        XMEN_OUT
        return;
    }
#ifndef XMEN_SIMPLE_BASE
    if (xmen_mip->base == 0)
    {
#endif
#ifndef XMEN_NOFULLWIDTH
        xg = x->w;
        yg = y->w;
        w0g = w0->w;
        if (x == y && xl > SQR_FASTER_THRESHOLD)
        /* extra hassle make it not    */
        /* worth it for small numbers */
        { /* fast squaring */
            for (i = 0; i < xl - 1; i++)
            { /* long multiplication */
                /* inline - substitutes for loop below */
                carry = 0;
                for (j = i + 1; j < xl; j++)
                { /* Only do above the diagonal */
                    dble.d = (XMen_LargeType)x->w[i] * x->w[j] + carry + w0->w[i + j];
                    w0->w[i + j] = dble.h[BYTE_BOT];
                    carry = dble.h[BYTE_TOP];
                }
                w0->w[xl + i] = carry;
            }
            w0->len = xl + xl - 1;
            xmen_padd(w0, w0, w0); /* double it */
            carry = 0;
            for (i = 0; i < xl; i++)
            { /* add in squared elements */
                ti = i + i;
                dble.d = (XMen_LargeType)x->w[i] * x->w[i] + carry + w0->w[ti];
                w0->w[ti] = dble.h[BYTE_BOT];
                carry = dble.h[BYTE_TOP];

                w0->w[ti + 1] += carry;
                if (w0->w[ti + 1] < carry)
                    carry = 1;
                else
                    carry = 0;
            }
        }
        else
            for (i = 0; i < xl; i++)
            { /* long multiplication */
                /* inline - substitutes for loop below */
                carry = 0;
                for (j = 0; j < yl; j++)
                { /* multiply each digit of y by x[i] */
                    dble.d = (XMen_LargeType)x->w[i] * y->w[j] + carry + w0->w[i + j];
                    w0->w[i + j] = dble.h[BYTE_BOT];
                    carry = dble.h[BYTE_TOP];
                }
                w0->w[yl + i] = carry;
            }
#endif
#ifndef XMEN_SIMPLE_BASE
    }
    else
    {
        if (x == y && xl > SQR_FASTER_THRESHOLD)
        { /* squaring can be done nearly twice as fast */
            for (i = 0; i < xl - 1; i++)
            { /* long multiplication */
                carry = 0;
                for (j = i + 1; j < xl; j++)
                { /* Only do above the diagonal */
                    dbled = (XMen_LargeType)x->w[i] * x->w[j] + w0->w[i + j] + carry;
#ifdef XMEN_FP_ROUNDING
                    carry = (XMen_SmallType)XMEN_LROUND(dbled * xmen_mip->inverse_base);
#else
                    if (xmen_mip->base == xmen_mip->base2)
                        carry = (XMen_SmallType)(dbled >> xmen_mip->lg2b);
                    else
                        carry = (XMen_SmallType)XMEN_LROUND(dbled / xmen_mip->base);
#endif
                    w0->w[i + j] = (XMen_SmallType)(dbled - (XMen_LargeType)carry * xmen_mip->base);
                }
                w0->w[xl + i] = carry;
            }
            w0->len = xl + xl - 1;
            xmen_padd(w0, w0, w0); /* double it */
            carry = 0;
            for (i = 0; i < xl; i++)
            { /* add in squared elements */
                ti = i + i;
                dbled = (XMen_LargeType)x->w[i] * x->w[i] + w0->w[ti] + carry;
#ifdef XMEN_FP_ROUNDING
                carry = (XMen_SmallType)XMEN_LROUND(dbled * xmen_mip->inverse_base);
#else
                if (xmen_mip->base == xmen_mip->base2)
                    carry = (XMen_SmallType)(dbled >> xmen_mip->lg2b);
                else
                    carry = (XMen_SmallType)XMEN_LROUND(dbled / xmen_mip->base);
#endif
                w0->w[ti] = (XMen_SmallType)(dbled - (XMen_LargeType)carry * xmen_mip->base);

                w0->w[ti + 1] += carry;
                carry = 0;
                if (w0->w[ti + 1] >= xmen_mip->base)
                {
                    carry = 1;
                    w0->w[ti + 1] -= xmen_mip->base;
                }
            }
        }
        else
            for (i = 0; i < xl; i++)
            { /* long multiplication */
                carry = 0;
                for (j = 0; j < yl; j++)
                { /* multiply each digit of y by x[i] */
                    dbled = (XMen_LargeType)x->w[i] * y->w[j] + w0->w[i + j] + carry;

#ifdef XMEN_FP_ROUNDING
                    carry = (XMen_SmallType)XMEN_LROUND(dbled * xmen_mip->inverse_base);
#else
                    if (xmen_mip->base == xmen_mip->base2)
                        carry = (XMen_SmallType)(dbled >> xmen_mip->lg2b);
                    else
                        carry = (XMen_SmallType)XMEN_LROUND(dbled / xmen_mip->base);
#endif
                    w0->w[i + j] = (XMen_SmallType)(dbled - (XMen_LargeType)carry * xmen_mip->base);
                }
                w0->w[yl + i] = carry;
            }
    }
#endif
    w0->len = (sz | (xl + yl)); /* set length and sign of result */

    xmen_lzero(w0);
    copy(w0, z);
    XMEN_OUT
}

/**
 * @description Initializes or sets a point on the active elliptic curve.
 *        If `x` or `y` is NULL, the function sets the point to the point at infinity.
 *        If `x == y`, it calculates a `y` coordinate if possible, using `cb` to suggest the least significant bit (LSB) of `y`.
 *        Otherwise, it checks the validity of the given `(x, y)` point while ignoring `cb`.
 *        The function returns TRUE if the point is valid, otherwise FALSE.
 * @param[in] x  : Pointer to the multi-precision integer representing the x-coordinate of the point. Can be NULL.
 * @param[in] y  : Pointer to the multi-precision integer representing the y-coordinate of the point. Can be NULL.
 * @param[in] cb : An integer suggesting the LSB of `y` when `x == y`. Ignored if `x != y`.
 * @param[out] p : Pointer to the `XMen_ECurvePointT` structure where the initialized point will be stored.
 * @return BOOL  : Returns TRUE if the point is successfully set and valid, otherwise FALSE.
 * @retval TRUE  : The point has been successfully set and is valid.
 * @retval FALSE : An error occurred, or the provided point is not valid on the elliptic curve.
 */
BOOL XMen_EpointSet(XMen_BigTPtr x, XMen_BigTPtr y, int cb, XMen_ECurvePointT *p)
{ /* initialise a point on active ecurve            *
   * if x or y == NULL, set to point at infinity    *
   * if x==y, a y co-ordinate is calculated - if    *
   * possible - and cb suggests LSB 0/1  of y       *
   * (which "decompresses" y). Otherwise, check     *
   * validity of given (x,y) point, ignoring cb.    *
   * Returns TRUE for valid point, otherwise FALSE. */

    BOOL valid;
    if (xmen_mip->ERNUM)
        return FALSE;

    XMEN_IN(97)

    if (x == NULL || y == NULL)
    {
        copy(xmen_mip->one, p->X);
        copy(xmen_mip->one, p->Y);
        p->marker = XMEN_EPOINT_INFINITY;
        XMEN_OUT
        return TRUE;
    }

    /* find x^3+Ax+B */

    nres(x, p->X);

    epoint_getrhs(p->X, xmen_mip->w3);

    valid = FALSE;

    if (x != y)
    { /* compare with y^2 */
        nres(y, p->Y);
        nres_modmult(p->Y, p->Y, xmen_mip->w1);

        if (XMen_Compare(xmen_mip->w1, xmen_mip->w3) == 0)
            valid = TRUE;
    }
    else
    { /* no y supplied - calculate one. Find square root */

        valid = nres_sqroot(xmen_mip->w3, p->Y);
        /* check LSB - have we got the right root? */
        redc(p->Y, xmen_mip->w1);
        if (remain(xmen_mip->w1, 2) != cb)
            xmen_psub(xmen_mip->modulus, p->Y, p->Y);
    }
    if (valid)
    {
        p->marker = XMEN_EPOINT_NORMALIZED;
        XMEN_OUT
        return TRUE;
    }

    XMEN_OUT
    return FALSE;
}

/**
 * @description Initializes the random number generation system.
 *        This function sets up the internal state of the random number generator using a provided seed value.
 *        It fills an initialization vector with values derived from the seed and then runs the generator multiple times
 *        to "warm-up" and stir the state, ensuring better randomness for subsequent calls.
 * @param[in] seed : An unsigned integer used as the seed to initialize the random number generator.
 * @return void
 */
void XMen_Irand(unsigned int seed)
{ /* initialise random number system */
    int i, in;
    unsigned int t, m = 1L;
    xmen_mip->borrow = 0L;
    xmen_mip->rndptr = 0;
    xmen_mip->ira[0] = seed;
    for (i = 1; i < NK; i++)
    { /* fill initialisation vector */
        in = (NV * i) % NK;
        xmen_mip->ira[in] = m;
        t = m;
        m = seed - m;
        seed = t;
    }
    for (i = 0; i < 1000; i++)
        brand(); /* "warm-up" & stir the generator */
}

/**
 * @description Generates a large random number within the specified range.
 *        This function creates a multi-precision integer `x` that is a random number satisfying `0 <= x < w`.
 *        It builds the random number piece by piece until it falls within the desired range, ensuring uniform distribution.
 * @param[in] w  : Pointer to the upper bound (exclusive) for the generated random number. Must be greater than zero.
 * @param[out] x : Pointer to the multi-precision integer where the generated random number will be stored.
 * @return void
 */
void XMen_BigRand(XMen_BigTPtr w, XMen_BigTPtr x)
{ /*  generate a big random number 0<=x<w  */
    int m;
    XMen_SmallType r;
    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(20)

    /*   decr(w,2,w);  */
    m = 0;
    zero(xmen_mip->w0);

    do
    { /* create big rand piece by piece */
        m++;
        xmen_mip->w0->len = m;
        r = brand();
        if (xmen_mip->base == 0)
            xmen_mip->w0->w[m - 1] = r;
        else
            xmen_mip->w0->w[m - 1] = XMEN_REMAIN(r, xmen_mip->base);
    } while (XMen_Compare(xmen_mip->w0, w) < 0);
    xmen_lzero(xmen_mip->w0);
    XMen_Divide(xmen_mip->w0, w, w);

    copy(xmen_mip->w0, x);
    XMEN_OUT
}

/**
 * @description Doubles a point on the active elliptic curve.
 *        This function computes the result of adding a point `p` to itself on the currently active elliptic curve,
 *        updating `p` with the resulting point. The operation is performed using either affine or projective coordinates,
 *        depending on the current setting, and handles special cases such as doubling the point at infinity or when the y-coordinate is zero.
 * @param[in,out] p : Pointer to the `XMen_ECurvePointT` structure representing the point to be doubled. The function updates this point with the result.
 * @return void
 */
void ecurve_double(XMen_ECurvePointT *p)
{ /* double XMen_ECurvePointT on active ecurve */
    if (xmen_mip->ERNUM)
        return;

    if (p->marker == XMEN_EPOINT_INFINITY)
    { /* 2 times infinity == infinity ! */
        return;
    }

    if (xmen_mip->coord == XMEN_AFFINE)
    { /* 2 sqrs, 1 mul, 1 div */
        if (EasyXMen_sm2_size(p->Y) == 0)
        { /* set to point at infinity */
            XMen_EpointSet(NULL, NULL, 0, p);
            return;
        }

        nres_modmult(p->X, p->X, xmen_mip->w8);      /* w8=x^2   */
        nres_premult(xmen_mip->w8, 3, xmen_mip->w8); /* w8=3*x^2 */
        if (xmen_abs(xmen_mip->Asize) == XMEN_TOOBIG)
            nres_modadd(xmen_mip->w8, xmen_mip->A, xmen_mip->w8);
        else
        {
            convert(xmen_mip->Asize, xmen_mip->w2);
            nres(xmen_mip->w2, xmen_mip->w2);
            nres_modadd(xmen_mip->w8, xmen_mip->w2, xmen_mip->w8);
        } /* w8=3*x^2+A */
        nres_premult(p->Y, 2, xmen_mip->w6); /* w6=2y */
        if (nres_moddiv(xmen_mip->w8, xmen_mip->w6, xmen_mip->w8) > 1)
        {
            XMen_EpointSet(NULL, NULL, 0, p);
            EasyXMen_berror(XMEN_ERR_COMPOSITE_MODULUS);
            return;
        }

        /* w8 is slope m on exit */

        nres_modmult(xmen_mip->w8, xmen_mip->w8, xmen_mip->w2); /* w2=m^2 */
        nres_premult(p->X, 2, xmen_mip->w1);
        nres_modsub(xmen_mip->w2, xmen_mip->w1, xmen_mip->w1); /* w1=m^2-2x */

        nres_modsub(p->X, xmen_mip->w1, xmen_mip->w2);
        nres_modmult(xmen_mip->w2, xmen_mip->w8, xmen_mip->w2);
        nres_modsub(xmen_mip->w2, p->Y, p->Y);
        copy(xmen_mip->w1, p->X);

        return;
    }

    if (EasyXMen_sm2_size(p->Y) == 0)
    { /* set to point at infinity */
        XMen_EpointSet(NULL, NULL, 0, p);
        return;
    }

    convert(1, xmen_mip->w1);
    if (xmen_abs(xmen_mip->Asize) < XMEN_TOOBIG)
    {
        if (xmen_mip->Asize != 0)
        {
            if (p->marker == XMEN_EPOINT_NORMALIZED)
                nres(xmen_mip->w1, xmen_mip->w6);
            else
                nres_modmult(p->Z, p->Z, xmen_mip->w6);
        }
        if (xmen_mip->Asize == (-3))
        { /* a is -3. Goody. 4 sqrs, 4 muls */
            nres_modsub(p->X, xmen_mip->w6, xmen_mip->w3);
            nres_modadd(p->X, xmen_mip->w6, xmen_mip->w8);
            nres_modmult(xmen_mip->w3, xmen_mip->w8, xmen_mip->w3);
            nres_modadd(xmen_mip->w3, xmen_mip->w3, xmen_mip->w8);
            nres_modadd(xmen_mip->w8, xmen_mip->w3, xmen_mip->w8);
        }
        else
        { /* a is small */
            if (xmen_mip->Asize != 0)
            { /* a is non zero! */
                nres_modmult(xmen_mip->w6, xmen_mip->w6, xmen_mip->w3);
                nres_premult(xmen_mip->w3, xmen_mip->Asize, xmen_mip->w3);
            }
            nres_modmult(p->X, p->X, xmen_mip->w1);
            nres_modadd(xmen_mip->w1, xmen_mip->w1, xmen_mip->w8);
            nres_modadd(xmen_mip->w8, xmen_mip->w1, xmen_mip->w8);
            if (xmen_mip->Asize != 0)
                nres_modadd(xmen_mip->w8, xmen_mip->w3, xmen_mip->w8);
        }
    }
    else
    { /* a is not special */
        if (p->marker == XMEN_EPOINT_NORMALIZED)
            nres(xmen_mip->w1, xmen_mip->w6);
        else
            nres_modmult(p->Z, p->Z, xmen_mip->w6);

        nres_modmult(xmen_mip->w6, xmen_mip->w6, xmen_mip->w3);
        nres_modmult(xmen_mip->w3, xmen_mip->A, xmen_mip->w3);
        nres_modmult(p->X, p->X, xmen_mip->w1);
        nres_modadd(xmen_mip->w1, xmen_mip->w1, xmen_mip->w8);
        nres_modadd(xmen_mip->w8, xmen_mip->w1, xmen_mip->w8);
        nres_modadd(xmen_mip->w8, xmen_mip->w3, xmen_mip->w8);
    }

    /* w8 contains numerator of slope 3x^2+A.z^4  *
     * denominator is now placed in Z             */

    nres_modmult(p->Y, p->Y, xmen_mip->w2);
    nres_modmult(p->X, xmen_mip->w2, xmen_mip->w3);
    nres_modadd(xmen_mip->w3, xmen_mip->w3, xmen_mip->w3);
    nres_modadd(xmen_mip->w3, xmen_mip->w3, xmen_mip->w3);
    nres_modmult(xmen_mip->w8, xmen_mip->w8, p->X);
    nres_modsub(p->X, xmen_mip->w3, p->X);
    nres_modsub(p->X, xmen_mip->w3, p->X);

    if (p->marker == XMEN_EPOINT_NORMALIZED)
        copy(p->Y, p->Z);
    else
        nres_modmult(p->Z, p->Y, p->Z);
    nres_modadd(p->Z, p->Z, p->Z);

    nres_modadd(xmen_mip->w2, xmen_mip->w2, xmen_mip->w7);
    nres_modmult(xmen_mip->w7, xmen_mip->w7, xmen_mip->w2);
    nres_modadd(xmen_mip->w2, xmen_mip->w2, xmen_mip->w2);
    nres_modsub(xmen_mip->w3, p->X, xmen_mip->w3);
    nres_modmult(xmen_mip->w8, xmen_mip->w3, p->Y);
    nres_modsub(p->Y, xmen_mip->w2, p->Y);

    p->marker = XMEN_EPOINT_GENERAL;
    return;
}

/**
 * @description Multiplies a point on the active elliptic curve by a scalar.
 *        This function calculates `pt = e * pa`, where `pa` is a point on the elliptic curve,
 *        and `e` is a multi-precision integer representing the scalar. The result is stored in `pt`.
 *        It uses precomputation tables for efficient computation and handles special cases, such as when `e` is zero or negative.
 * @param[in] e   : Pointer to the multi-precision integer representing the scalar value to multiply by.
 * @param[in] pa  : Pointer to the `XMen_ECurvePointT` structure representing the base point on the elliptic curve.
 * @param[out] pt : Pointer to the `XMen_ECurvePointT` structure where the resulting point will be stored.
 * @return int    : Returns the number of additions performed during the multiplication process.
 * @retval >=0    : Number of additions performed.
 */
int XMen_EcurveMult(XMen_BigTPtr e, XMen_ECurvePointT *pa, XMen_ECurvePointT *pt)
{ /* pt=e*pa; */
    int i, j, n, nb, nbs, nzs, nadds;
    XMen_ECurvePointT *table[XMEN_ECC_STORE_N];
    XMen_BigTPtr work[XMEN_ECC_STORE_N];

    char *mem;
    char *mem1;

#ifndef XMEN_ALWAYS_BINARY
    XMen_ECurvePointT *p;
    int ce, ch;
#endif
    if (xmen_mip->ERNUM)
        return 0;

    XMEN_IN(95)
    if (EasyXMen_sm2_size(e) == 0)
    { /* multiplied by 0 */
        XMen_EpointSet(NULL, NULL, 0, pt);
        XMEN_OUT
        return 0;
    }
    copy(e, xmen_mip->w9);
    /*    epoint_norm(pa); */
    epoint_copy(pa, pt);

    if (EasyXMen_sm2_size(xmen_mip->w9) < 0)
    { /* pt = -pt */
        negify(xmen_mip->w9, xmen_mip->w9);
        epoint_negate(pt);
    }

    if (EasyXMen_sm2_size(xmen_mip->w9) == 1)
    {
        XMEN_OUT
        return 0;
    }

    premult(xmen_mip->w9, 3, xmen_mip->w10); /* h=3*e */

#ifndef XMEN_STATIC
#ifndef XMEN_ALWAYS_BINARY
    if (xmen_mip->base == xmen_mip->base2)
    {
#endif
#endif

        mem = (char *)ecp_memalloc(XMEN_ECC_STORE_N);
        mem1 = (char *)EasyXMen_memalloc(XMEN_ECC_STORE_N);

        for (i = 0; i <= XMEN_ECC_STORE_N - 1; i++)
        {
            table[i] = epoint_init_mem(mem, i);
            work[i] = mirvar_mem(mem1, i);
        }

        epoint_copy(pt, table[0]);
        epoint_copy(table[0], table[XMEN_ECC_STORE_N - 1]);
        ecurve_double(table[XMEN_ECC_STORE_N - 1]);
        /*   epoint_norm(table[XMEN_ECC_STORE_N-1]); */

        for (i = 1; i < XMEN_ECC_STORE_N - 1; i++)
        { /* precomputation */
            epoint_copy(table[i - 1], table[i]);
            ecurve_add(table[XMEN_ECC_STORE_N - 1], table[i]);
        }
        ecurve_add(table[XMEN_ECC_STORE_N - 2], table[XMEN_ECC_STORE_N - 1]);

        epoint_multi_norm(XMEN_ECC_STORE_N, work, table); // madezheng

        nb = logb2(xmen_mip->w10);
        nadds = 0;
        XMen_EpointSet(NULL, NULL, 0, pt);
        for (i = nb - 1; i >= 1;)
        { /* add/subtract */
            if (xmen_mip->user != NULL)
                (*xmen_mip->user)();
            n = xmen_naf_window(xmen_mip->w9, xmen_mip->w10, i, &nbs, &nzs, XMEN_ECC_STORE_N);
            for (j = 0; j < nbs; j++)
                ecurve_double(pt);
            if (n > 0)
            {
                ecurve_add(table[n / 2], pt);
                nadds++;
            }
            if (n < 0)
            {
                ecurve_sub(table[(-n) / 2], pt);
                nadds++;
            }
            i -= nbs;
            if (nzs)
            {
                for (j = 0; j < nzs; j++)
                    ecurve_double(pt);
                i -= nzs;
            }
        }

        ecp_memkill(mem, XMEN_ECC_STORE_N);
        memkill(mem1, XMEN_ECC_STORE_N);

#ifndef XMEN_STATIC
#ifndef XMEN_ALWAYS_BINARY
    }
    else
    {
        mem = (char *)ecp_memalloc(1);
        p = epoint_init_mem(mem, 0);
        epoint_norm(pt);
        epoint_copy(pt, p);

        nadds = 0;
        expb2(logb2(xmen_mip->w10) - 1, xmen_mip->w11);
        xmen_psub(xmen_mip->w10, xmen_mip->w11, xmen_mip->w10);
        subdiv(xmen_mip->w11, 2, xmen_mip->w11);
        while (EasyXMen_sm2_size(xmen_mip->w11) > 1)
        { /* add/subtract method */
            if (xmen_mip->user != NULL)
                (*xmen_mip->user)();

            ecurve_double(pt);
            ce = XMen_Compare(xmen_mip->w9, xmen_mip->w11);  /* e(i)=1? */
            ch = XMen_Compare(xmen_mip->w10, xmen_mip->w11); /* h(i)=1? */
            if (ch >= 0)
            { /* h(i)=1 */
                if (ce < 0)
                {
                    ecurve_add(p, pt);
                    nadds++;
                }
                xmen_psub(xmen_mip->w10, xmen_mip->w11, xmen_mip->w10);
            }
            if (ce >= 0)
            { /* e(i)=1 */
                if (ch < 0)
                {
                    ecurve_sub(p, pt);
                    nadds++;
                }
                xmen_psub(xmen_mip->w9, xmen_mip->w11, xmen_mip->w9);
            }
            subdiv(xmen_mip->w11, 2, xmen_mip->w11);
        }
        ecp_memkill(mem, 1);
    }
#endif
#endif
    XMEN_OUT
    return nadds;
}

/**
 * @description Retrieves the coordinates of an elliptic curve point in affine, normalized form.
 *        This function converts the coordinates of a point `p` from projective or Montgomery form to affine form,
 *        and stores them in multi-precision integers `x` and `y`. If `x == y`, it supplies only the `x` coordinate.
 *        The return value is the least significant bit (LSB) of `y`, which can be useful for point compression.
 * @param[in] p  : Pointer to the `XMen_ECurvePointT` structure representing the point on the elliptic curve.
 * @param[out] x : Pointer to the multi-precision integer where the x-coordinate will be stored.
 * @param[out] y : Pointer to the multi-precision integer where the y-coordinate will be stored. If `x == y`, only `x` is set.
 * @return int   : Returns the least significant bit (LSB) of the y-coordinate.
 * @retval 0     : The point is at infinity, both `x` and `y` are set to zero.
 * @retval -1    : Normalization of the point was not possible.
 * @retval >=0   : The least significant bit (LSB) of the y-coordinate.
 */
int XMen_EpointGet(XMen_ECurvePointT *p, XMen_BigTPtr x, XMen_BigTPtr y)
{ /* Get point co-ordinates in affine, normal form       *
   * (converted from projective, Montgomery form)        *
   * if x==y, supplies x only. Return value is Least     *
   * Significant Bit of y (useful for point compression) */

    int lsb;

    if (p->marker == XMEN_EPOINT_INFINITY)
    {
        zero(x);
        zero(y);
        return 0;
    }
    if (xmen_mip->ERNUM)
        return 0;

    XMEN_IN(98)

    if (!epoint_norm(p))
    { /* not possible ! */
        XMEN_OUT
        return (-1);
    }

    redc(p->X, x);
    redc(p->Y, xmen_mip->w1);

    if (x != y)
        copy(xmen_mip->w1, y);
    lsb = remain(xmen_mip->w1, 2);
    XMEN_OUT
    return lsb;
}

/* @note:convert big into octet string */
/**
 * @description Converts a multi-precision integer to a byte array.
 *        This function converts the value of a multi-precision integer `x` into a byte array,
 *        storing it in the memory pointed to by `ptr`. It supports both justified and non-justified output formats.
 *        The function ensures that the byte array does not exceed the specified maximum length `max`.
 * @param[in] max       : The maximum number of bytes to write to `ptr`. If `justify` is TRUE, the output will be padded with leading zeros to this length.
 * @param[in] x         : Pointer to the multi-precision integer to convert.
 * @param[out] ptr      : Pointer to the buffer where the resulting byte array will be stored.
 * @param[in] justify   : A boolean flag indicating whether the output should be right-justified within `max` bytes.
 * @return int          : Returns the number of bytes written to `ptr`.
 * @retval 0            : Conversion failed due to an error (e.g., input size is zero or invalid parameters).
 * @retval >0           : The number of bytes written to `ptr`.
 */
int XMen_Big2Bytes(int max, XMen_BigTPtr x, char *ptr, BOOL justify)
{
    int i, j, r, m, n, len, start;
    unsigned int dig;
    unsigned char ch;
    XMen_SmallType wrd;

    if (xmen_mip->ERNUM)
        return 0;

    if (EasyXMen_sm2_size(x) == 0 || (max <= 0 && justify))
        return 0;
    XMEN_IN(141);

    xmen_lzero(x); /* should not be needed.... */
#ifndef XMEN_SIMPLE_BASE
    if (xmen_mip->base == 0)
    {
#endif
#ifndef XMEN_NOFULLWIDTH
        m = MIRACL / 8;
        n = (int)(x->len & XMEN_OBITS);
        n--;
        len = n * m;
        wrd = x->w[n]; /* most significant */
        r = 0;
        while (wrd != (XMen_SmallType)0)
        {
            r++;
            wrd >>= 8;
            len++;
        }
        r %= m;

        if (max > 0 && len > max)
        {
            EasyXMen_berror(XMEN_ERR_TOO_BIG);
            XMEN_OUT
            return 0;
        }

        if (justify)
        {
            start = max - len;
            for (i = 0; i < start; i++)
                ptr[i] = 0;
        }
        else
            start = 0;

        if (r != 0)
        {
            wrd = x->w[n--];
            for (i = r - 1; i >= 0; i--)
            {
                ptr[start + i] = (char)(wrd & 0xFF);
                wrd >>= 8;
            }
        }

        for (i = r; i < len; i += m)
        {
            wrd = x->w[n--];
            for (j = m - 1; j >= 0; j--)
            {
                ptr[start + i + j] = (char)(wrd & 0xFF);
                wrd >>= 8;
            }
        }
#endif
#ifndef XMEN_SIMPLE_BASE
    }
    else
    {
        copy(x, xmen_mip->w1);
        for (len = 0;; len++)
        {
            if (xmen_mip->ERNUM)
                break;

            if (EasyXMen_sm2_size(xmen_mip->w1) == 0)
            {
                if (justify)
                {
                    if (len == max)
                        break;
                }
                else
                    break;
            }

            if (max > 0 && len >= max)
            {
                EasyXMen_berror(XMEN_ERR_TOO_BIG);
                XMEN_OUT
                return 0;
            }
            dig = (unsigned int)subdiv(xmen_mip->w1, 256, xmen_mip->w1);
            ch = XMEN_TOBYTE(dig);
            for (i = len; i > 0; i--)
                ptr[i] = ptr[i - 1];
            ptr[0] = XMEN_TOBYTE(ch);
        }
    }
#endif
    XMEN_OUT
    if (justify)
        return max;
    else
        return len;
}

#ifndef XMEN_STATIC
/**
 * @description Destroys a multi-precision integer variable, setting it to zero and freeing its allocated memory.
 *        This function ensures that any sensitive data stored in the multi-precision integer `x` is securely erased,
 *        and that the memory used by `x` is deallocated.
 * @param[in,out] x : Pointer to the multi-precision integer to be destroyed. If `x` is NULL, the function does nothing.
 * @return void
 */
void XMen_Mirkill(XMen_BigTPtr x)
{ /* kill a big/flash variable, that is set it to zero
     and free its memory */
    if (x == NULL)
        return;
    zero(x);
    EasyXMen_free(x);
}

#endif
/**
 * @description Cleans up and frees an elliptic curve point.
 *        This function securely erases the coordinates of the elliptic curve point `p` by setting them to zero,
 *        and then deallocates the memory used by `p`. If `p` is NULL, the function does nothing.
 * @param[in,out] p : Pointer to the `XMen_ECurvePointT` structure representing the elliptic curve point to be freed.
 *                    The function sets the coordinates of `p` to zero and frees the memory allocated for `p`.
 * @return void
 */
void XMen_EpointFree(XMen_ECurvePointT *p)
{ /* clean up point */

    if (p == NULL)
        return;
    zero(p->X);
    zero(p->Y);
    if (p->marker == XMEN_EPOINT_GENERAL)
        zero(p->Z);
    EasyXMen_free(p);
}

/**
 * @description Cleans up the environment after using multi-precision integer and elliptic curve operations.
 *        This function performs necessary cleanup actions, such as deactivating the current context,
 *        freeing allocated memory, resetting error conditions, and clearing any sensitive data.
 *        It ensures that resources are properly released and that no residual data remains in memory.
 * @return void
 */
void XMen_Mirexit(void)
{ /* clean up after XMen_SM2Type */

    int i;
    xmen_mip->ERCON = FALSE;
    xmen_mip->active = OFF;
    memkill(xmen_mip->workspace, XMEN_SPACES);
#ifndef XMEN_NO_RAND
    for (i = 0; i < NK; i++)
        xmen_mip->ira[i] = 0L;
#endif
#ifndef XMEN_STATIC
#ifndef XMEN_SIMPLE_IO
    set_io_buffer_size(0);
#endif
    if (xmen_mip->PRIMES != NULL)
        EasyXMen_free(xmen_mip->PRIMES);
#else
#endif

#ifndef XMEN_STATIC
    EasyXMen_free(xmen_mip);
#endif

#ifndef XMEN_STATIC
    xmen_mip = NULL;
#endif
}

/**
 * @description Computes a linear combination of two elliptic curve points.
 *        This function calculates `pt = e*p + ea*pa`, where `p` and `pa` are points on the elliptic curve,
 *        and `e` and `ea` are multi-precision integers representing the scalar values to multiply by each point.
 *        The result is stored in `pt`.
 * @param[in] e   : Pointer to the multi-precision integer representing the scalar value to multiply with point `p`.
 * @param[in] p   : Pointer to the `XMen_ECurvePointT` structure representing the first point on the elliptic curve.
 * @param[in] ea  : Pointer to the multi-precision integer representing the scalar value to multiply with point `pa`.
 * @param[in] pa  : Pointer to the `XMen_ECurvePointT` structure representing the second point on the elliptic curve.
 * @param[out] pt : Pointer to the `XMen_ECurvePointT` structure where the resulting point will be stored.
 * @return void
 */
void XMen_EcurveMult2(XMen_BigTPtr e, XMen_ECurvePointT *p, XMen_BigTPtr ea, XMen_ECurvePointT *pa, XMen_ECurvePointT *pt)
{ /* pt=e*p+ea*pa; */
    int e1, h1, e2, h2, bb;
    XMen_ECurvePointT *p1, *p2, *ps[2];
    char *mem;

    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(103)

    if (EasyXMen_sm2_size(e) == 0)
    {
        XMen_EcurveMult(ea, pa, pt);
        XMEN_OUT
        return;
    }
    mem = (char *)ecp_memalloc(4);
    p2 = epoint_init_mem(mem, 0);
    p1 = epoint_init_mem(mem, 1);
    ps[0] = epoint_init_mem(mem, 2);
    ps[1] = epoint_init_mem(mem, 3);

    epoint_norm(pa);
    epoint_copy(pa, p2);
    copy(ea, xmen_mip->w9);
    if (EasyXMen_sm2_size(xmen_mip->w9) < 0)
    { /* p2 = -p2 */
        negify(xmen_mip->w9, xmen_mip->w9);
        epoint_negate(p2);
    }

    epoint_norm(p);
    epoint_copy(p, p1);
    copy(e, xmen_mip->w12);
    if (EasyXMen_sm2_size(xmen_mip->w12) < 0)
    { /* p1= -p1 */
        negify(xmen_mip->w12, xmen_mip->w12);
        epoint_negate(p1);
    }

    XMen_EpointSet(NULL, NULL, 0, pt);    /* pt=0 */
    ecurve_add_sub(p1, p2, ps[0], ps[1]); /* only one inversion! ps[0]=p1+p2, ps[1]=p1-p2 */

    xmen_jsf(xmen_mip->w9, xmen_mip->w12, xmen_mip->w10, xmen_mip->w9, xmen_mip->w13, xmen_mip->w12);

    /*    To use a simple NAF instead, substitute this for the JSF
            premult(xmen_mip->w9,3,xmen_mip->w10);      3*ea
            premult(xmen_mip->w12,3,xmen_mip->w13);     3*e
    */

#ifndef XMEN_ALWAYS_BINARY
    if (xmen_mip->base == xmen_mip->base2)
    {
#endif
        if (XMen_Compare(xmen_mip->w10, xmen_mip->w13) >= 0)
            bb = logb2(xmen_mip->w10) - 1;
        else
            bb = logb2(xmen_mip->w13) - 1;

        while (bb >= 0) /* for the simple NAF, this should be 1 */
        {
            if (xmen_mip->user != NULL)
                (*xmen_mip->user)();
            ecurve_double(pt);

            e1 = h1 = e2 = h2 = 0;
            if (xmen_testbit(xmen_mip->w9, bb))
                e2 = 1;
            if (xmen_testbit(xmen_mip->w10, bb))
                h2 = 1;
            if (xmen_testbit(xmen_mip->w12, bb))
                e1 = 1;
            if (xmen_testbit(xmen_mip->w13, bb))
                h1 = 1;

            if (e1 != h1)
            {
                if (e2 == h2)
                {
                    if (h1 == 1)
                        ecurve_add(p1, pt);
                    else
                        ecurve_sub(p1, pt);
                }
                else
                {
                    if (h1 == 1)
                    {
                        if (h2 == 1)
                            ecurve_add(ps[0], pt);
                        else
                            ecurve_add(ps[1], pt);
                    }
                    else
                    {
                        if (h2 == 1)
                            ecurve_sub(ps[1], pt);
                        else
                            ecurve_sub(ps[0], pt);
                    }
                }
            }
            else if (e2 != h2)
            {
                if (h2 == 1)
                    ecurve_add(p2, pt);
                else
                    ecurve_sub(p2, pt);
            }
            bb -= 1;
        }
#ifndef XMEN_ALWAYS_BINARY
    }
    else
    {
        if (XMen_Compare(xmen_mip->w10, xmen_mip->w13) >= 0)
            expb2(logb2(xmen_mip->w10) - 1, xmen_mip->w11);
        else
            expb2(logb2(xmen_mip->w13) - 1, xmen_mip->w11);

        while (EasyXMen_sm2_size(xmen_mip->w11) > 0) /* for the NAF, this should be 1 */
        {                                            /* XMen_Add/XMen_Subtract method */
            if (xmen_mip->user != NULL)
                (*xmen_mip->user)();

            ecurve_double(pt);

            e1 = h1 = e2 = h2 = 0;
            if (XMen_Compare(xmen_mip->w9, xmen_mip->w11) >= 0)
            { /* e1(i)=1? */
                e2 = 1;
                xmen_psub(xmen_mip->w9, xmen_mip->w11, xmen_mip->w9);
            }
            if (XMen_Compare(xmen_mip->w10, xmen_mip->w11) >= 0)
            { /* h1(i)=1? */
                h2 = 1;
                xmen_psub(xmen_mip->w10, xmen_mip->w11, xmen_mip->w10);
            }
            if (XMen_Compare(xmen_mip->w12, xmen_mip->w11) >= 0)
            { /* e2(i)=1? */
                e1 = 1;
                xmen_psub(xmen_mip->w12, xmen_mip->w11, xmen_mip->w12);
            }
            if (XMen_Compare(xmen_mip->w13, xmen_mip->w11) >= 0)
            { /* h2(i)=1? */
                h1 = 1;
                xmen_psub(xmen_mip->w13, xmen_mip->w11, xmen_mip->w13);
            }

            if (e1 != h1)
            {
                if (e2 == h2)
                {
                    if (h1 == 1)
                        ecurve_add(p1, pt);
                    else
                        ecurve_sub(p1, pt);
                }
                else
                {
                    if (h1 == 1)
                    {
                        if (h2 == 1)
                            ecurve_add(ps[0], pt);
                        else
                            ecurve_add(ps[1], pt);
                    }
                    else
                    {
                        if (h2 == 1)
                            ecurve_sub(ps[1], pt);
                        else
                            ecurve_sub(ps[0], pt);
                    }
                }
            }
            else if (e2 != h2)
            {
                if (h2 == 1)
                    ecurve_add(p2, pt);
                else
                    ecurve_sub(p2, pt);
            }

            subdiv(xmen_mip->w11, 2, xmen_mip->w11);
        }
    }
#endif
    ecp_memkill(mem, 4);
    XMEN_OUT
}

/*------------------------------Internal function define.------------------------------*/
/**
 * @description Retrieves the current multi-precision integer processing context.
 *        This function returns a pointer to the `XMen_SM2Type` structure that holds the current context for multi-precision integer operations.
 * @return XMen_SM2Type* : Pointer to the current multi-precision integer processing context.
 */
static XMen_SM2Type *get_mip()
{
    return (XMen_SM2Type *)xmen_mip;
}

/**
 * @description Converts a signed integer to multi-precision integer format.
 *        This function takes a signed integer `n` and converts it into the big number format used by the library,
 *        storing the result in the multi-precision integer `x`. If `n` is negative, its absolute value is converted,
 *        and the sign is preserved in the resulting multi-precision integer.
 * @param[in] n  : The signed integer to convert.
 * @param[out] x : Pointer to the multi-precision integer where the converted value will be stored.
 * @return void
 */
static void convert(int n, XMen_BigTPtr x)
{ /*  convert signed integer n to big number format  */
    unsigned int s;

    if (n == 0)
    {
        zero(x);
        return;
    }
    s = 0;
    if (n < 0)
    {
        s = XMEN_MSBIT;
        n = (-n);
    }
    uconvert((unsigned int)n, x);
    x->len |= s;
}

/**
 * @description Allocates and aligns memory for a multi-precision integer variable from a pre-allocated memory block.
 *        This function sets up a multi-precision integer `x` within a specified memory area `mem`, ensuring proper alignment
 *        according to the system's requirements. The allocated variable is intended to be used as a flash-based multi-precision integer.
 * @param[in] mem : Pointer to the pre-allocated memory block where the multi-precision integer will be set up.
 * @param[in] index : Index into the memory block to determine the starting position of the new variable.
 * @param[in] sz : Size parameter that affects the allocation size and alignment of the multi-precision integer.
 * @return flash : A pointer to the newly allocated and aligned multi-precision integer variable.
 */
static flash mirvar_mem_variable(char *mem, int index, int sz)
{
    flash x;
    int align;
    char *ptr;
    int offset, r;

    /* alignment */
    offset = 0;
    r = (unsigned long)mem % XMEN_SL;
    if (r > 0)
        offset = XMEN_SL - r;

    x = (XMen_BigTPtr)&mem[offset + xmen_size(sz) * index];
    ptr = (char *)&x->w;
    align = (unsigned long)(ptr + sizeof(XMen_SmallType *)) % sizeof(XMen_SmallType);
    x->w = (XMen_SmallType *)(ptr + sizeof(XMen_SmallType *) + sizeof(XMen_SmallType) - align);

    return x;
}

/**
 * @description Initializes a multi-precision integer variable from pre-allocated memory.
 *        This function sets up a multi-precision integer `x` within a specified memory area `mem`, using the given `index`
 *        to determine its position. It relies on `mirvar_mem_variable` for the actual allocation and alignment of the variable.
 * @param[in] mem : Pointer to the pre-allocated memory block where the multi-precision integer will be initialized.
 * @param[in] index : Index into the memory block to determine the starting position of the new variable.
 * @return flash : A pointer to the newly initialized multi-precision integer variable, or NULL if an error occurs.
 * @retval NULL : An error occurred during initialization, possibly due to a previous error state.
 */
static flash mirvar_mem(char *mem, int index)
{ /* initialize big/flash number from pre-allocated memory */

    if (xmen_mip->ERNUM)
        return NULL;

    return mirvar_mem_variable(mem, index, xmen_mip->nib - 1);
}

/**
 * @description Determines the size of a multi-precision integer and attempts to convert it to an integer if possible.
 *        This function checks the length and value of the multi-precision integer `x` and returns its size as an integer.
 *        If `x` is negative, the returned value will also be negative. If `x` cannot be represented as an integer,
 *        a special value `XMEN_TOOBIG` is returned to indicate overflow.
 * @param[in] x : Pointer to the multi-precision integer whose size is to be determined.
 * @return int  : The size of the multi-precision integer as an integer value, or 0 if `x` is NULL or has zero length.
 *                Returns `-n` if `x` is negative, where `n` is the absolute value of `x`.
 *                Returns `XMEN_TOOBIG` if the value of `x` exceeds the representable range of an integer.
 * @retval 0    : If `x` is NULL or has zero length.
 */
static int EasyXMen_sm2_size(XMen_BigTPtr x)
{ /*  get size of big number;  convert to *
   *  integer - if possible               */
    int n, m;
    unsigned int s;
    if (x == NULL)
        return 0;
    s = (x->len & XMEN_MSBIT);
    m = (int)(x->len & XMEN_OBITS);
    if (m == 0)
        return 0;
    if (m == 1 && x->w[0] < (XMen_SmallType)XMEN_TOOBIG)
        n = (int)x->w[0];
    else
        n = XMEN_TOOBIG;
    if (s == XMEN_MSBIT)
        return (-n);
    return n;
}

/**
 * @description Multiplies an n-residue by a small ordinary integer.
 *        This function performs multiplication of the multi-precision integer `x` (interpreted as an n-residue) by a small integer `k`,
 *        storing the result in `w`. If `k` is negative, the result is negated. The function handles special cases for small values of `k`
 *        using optimized operations and falls back to a general-purpose multiplication algorithm for larger values.
 * @param[in] x   : Pointer to the multi-precision integer (n-residue) to be multiplied.
 * @param[in] k   : The small integer multiplier.
 * @param[out] w  : Pointer to the multi-precision integer where the result will be stored.
 * @return void
 */
static void nres_premult(XMen_BigTPtr x, int k, XMen_BigTPtr w)
{ /* multiply n-residue by small ordinary integer */
    int sign = 0;
    if (k == 0)
    {
        zero(w);
        return;
    }
    if (k < 0)
    {
        k = -k;
        sign = 1;
    }
    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(102)

    if (k <= 6)
    {
        switch (k)
        {
        case 1:
            copy(x, w);
            break;
        case 2:
            nres_modadd(x, x, w);
            break;
        case 3:
            nres_modadd(x, x, xmen_mip->w0);
            nres_modadd(x, xmen_mip->w0, w);
            break;
        case 4:
            nres_modadd(x, x, w);
            nres_modadd(w, w, w);
            break;
        case 5:
            nres_modadd(x, x, xmen_mip->w0);
            nres_modadd(xmen_mip->w0, xmen_mip->w0, xmen_mip->w0);
            nres_modadd(x, xmen_mip->w0, w);
            break;
        case 6:
            nres_modadd(x, x, w);
            nres_modadd(w, w, xmen_mip->w0);
            nres_modadd(w, xmen_mip->w0, w);
            break;
        }
        if (sign == 1)
            nres_negate(w, w);
        XMEN_OUT
        return;
    }

    xmen_pmul(x, (XMen_SmallType)k, xmen_mip->w0);
#ifdef XMEN_COMBA
#else
    XMen_Divide(xmen_mip->w0, xmen_mip->modulus, xmen_mip->modulus);
    copy(xmen_mip->w0, w);
#endif

    if (sign == 1)
        nres_negate(w, w);

    XMEN_OUT
}

/**
 * @description Performs modular addition of two multi-precision integers.
 *        This function adds the multi-precision integers `x` and `y` modulo a specified modulus,
 *        storing the result in `w`. It checks for an overflow condition and adjusts the result accordingly.
 *        If the COMBA optimization is active, it uses the optimized `comba_modadd` function for the operation.
 * @param[in] x   : Pointer to the first multi-precision integer operand.
 * @param[in] y   : Pointer to the second multi-precision integer operand.
 * @param[out] w  : Pointer to the multi-precision integer where the result will be stored.
 * @return void
 */
static void nres_modadd(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr w)
{ /* modular addition */
#ifdef XMEN_COMBA

    if (xmen_mip->ACTIVE)
    {
        comba_modadd(x, y, w);
        return;
    }
    else
    {
#endif
        if (xmen_mip->ERNUM)
            return;

        XMEN_IN(90)
        xmen_padd(x, y, w);
        if (XMen_Compare(w, xmen_mip->modulus) >= 0)
            xmen_psub(w, xmen_mip->modulus, w);

        XMEN_OUT
#ifdef XMEN_COMBA
    }
#endif
}

/**
 * @description Computes the Jacobi symbol (a/n) for positive odd n.
 *        This function calculates the Jacobi symbol of two multi-precision integers `a` and `n`, where `n` must be a positive odd integer.
 *        The Jacobi symbol is an extension of the Legendre symbol used in number theory, indicating whether `a` is a quadratic residue modulo `n`.
 * @param[in] a   : Pointer to the first multi-precision integer operand.
 * @param[in] n   : Pointer to the second multi-precision integer operand, which should be positive and odd.
 * @return int    : The computed Jacobi symbol value, or 0 if `n` is not positive and odd or if an error occurs.
 * @retval 1      : If `a` is a quadratic residue modulo `n`.
 * @retval -1     : If `a` is a non-quadratic residue modulo `n`.
 * @retval 0      : If `a` is divisible by `n`, or if `n` is not positive and odd, or if an error occurs.
 */
static int jack(XMen_BigTPtr a, XMen_BigTPtr n)
{ /* find jacobi symbol (a/n), for positive odd n */
    XMen_BigTPtr w;
    int nm8, onm8, t;
    if (xmen_mip->ERNUM || EasyXMen_sm2_size(a) == 0 || EasyXMen_sm2_size(n) < 1)
        return 0;
    XMEN_IN(3)

    t = 1;
    copy(n, xmen_mip->w2);
    nm8 = remain(xmen_mip->w2, 8);
    if (nm8 % 2 == 0)
    {
        XMEN_OUT
        return 0;
    }

    if (EasyXMen_sm2_size(a) < 0)
    {
        if (nm8 % 4 == 3)
            t = -1;
        negify(a, xmen_mip->w1);
    }
    else
        copy(a, xmen_mip->w1);

    while (EasyXMen_sm2_size(xmen_mip->w1) != 0)
    {
        while (remain(xmen_mip->w1, 2) == 0)
        {
            subdiv(xmen_mip->w1, 2, xmen_mip->w1);
            if (nm8 == 3 || nm8 == 5)
                t = -t;
        }
        if (XMen_Compare(xmen_mip->w1, xmen_mip->w2) < 0)
        {
            onm8 = nm8;
            w = xmen_mip->w1;
            xmen_mip->w1 = xmen_mip->w2;
            xmen_mip->w2 = w;
            nm8 = remain(xmen_mip->w2, 8);
            if (onm8 % 4 == 3 && nm8 % 4 == 3)
                t = -t;
        }
        xmen_psub(xmen_mip->w1, xmen_mip->w2, xmen_mip->w1);
        subdiv(xmen_mip->w1, 2, xmen_mip->w1);

        if (nm8 == 3 || nm8 == 5)
            t = -t;
    }

    XMEN_OUT
    if (EasyXMen_sm2_size(xmen_mip->w2) == 1)
        return t;
    return 0;
}

/**
 * @description Extracts the denominator of a multi-precision rational number.
 *        This function extracts the denominator part of the given flash number `x` and stores it in the multi-precision integer `y`.
 *        If `x` is an integer (i.e., its denominator is 1), then `y` is set to 1. Otherwise, the denominator is extracted from `x` and stored in `y`.
 * @param[in] x   : Pointer to the flash number from which to extract the denominator.
 * @param[out] y  : Pointer to the multi-precision integer where the extracted denominator will be stored.
 * @return void
 */
static void denom(flash x, XMen_BigTPtr y)
{ /* extract denominator of x */
    int i, ln, ld;
    unsigned int ly;

    if (xmen_mip->ERNUM)
        return;
    if (!xmen_notint(x))
    {
        convert(1, y);
        return;
    }
    ly = (x->len & XMEN_OBITS);
    ln = (int)(ly & XMEN_MSK);
    ld = (int)((ly >> XMEN_BTS) & XMEN_MSK);
    for (i = 0; i < ld; i++)
        y->w[i] = x->w[ln + i];
    if (x == y)
        for (i = 0; i < ln; i++)
            y->w[ld + i] = 0;
    else
        for (i = ld; i < xmen_lent(y); i++)
            y->w[i] = 0;
    y->len = ld;
}

/**
 * @description Asserts the sign of a multi-precision or flash number.
 *        This function sets the sign bit of the given flash number `x` based on the input parameter `s`.
 *        If `s` is negative, the sign bit is set to indicate a negative number; otherwise, it is cleared.
 * @param[in] s   : An integer indicating the desired sign of the number (-1 for negative, 0 or positive for non-negative).
 * @param[inout] x: Pointer to the flash number whose sign is to be asserted.
 * @return void
 */
static void insign(int s, flash x)
{ /* assert sign of big/flash number */
    if (x->len == 0)
        return;
    if (s < 0)
        x->len |= XMEN_MSBIT;
    else
        x->len &= XMEN_OBITS;
}

/**
 * @description Shifts the bits of a small multi-precision integer by a specified number of positions.
 *        This function shifts the bits of the `XMen_SmallType` value `x` left if `n` is positive, or right if `n` is negative,
 *        and returns the result. If `n` is zero, the original value of `x` is returned unchanged.
 * @param[in] x   : The `XMen_SmallType` value whose bits are to be shifted.
 * @param[in] n   : The number of positions to shift the bits of `x`. Positive values shift left, negative values shift right.
 * @return XMen_SmallType : The result of shifting the bits of `x` by `n` positions.
 */
static XMen_SmallType xmen_shiftbits(XMen_SmallType x, int n)
{
    if (n == 0)
        return x;
    if (n > 0)
        x <<= n;
    else
        x >>= (-n);
    return x;
}

/**
 * @description Creates a floating-slash number from a big integer numerator and denominator.
 *        This function constructs a flash number `x` that represents the rational value `n/d`, where `n` is the numerator and `d` is the denominator.
 *        It performs various checks to ensure the operation's validity, including checking for overflow conditions and ensuring the denominator is not zero.
 * @param[in] n   : Pointer to the multi-precision integer representing the numerator.
 * @param[in] d   : Pointer to the multi-precision integer representing the denominator.
 * @param[out] x  : Pointer to the flash number where the result will be stored.
 * @return void
 */
static void fpack(XMen_BigTPtr n, XMen_BigTPtr d, flash x)
{ /* create floating-slash number x=n/d from *
   * big integer numerator and denominator   */
    unsigned int s;
    int i, ld, ln;

    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(31)

    ld = (int)(d->len & XMEN_OBITS);
    if (ld == 0)
        EasyXMen_berror(XMEN_ERR_FLASH_OVERFLOW);
    if (ld == 1 && d->w[0] == 1)
        ld = 0;
    if (x == d)
        EasyXMen_berror(XMEN_ERR_BAD_PARAMETERS);
    if (xmen_notint(n) || xmen_notint(d))
        EasyXMen_berror(XMEN_ERR_INT_OP);
    s = (n->len & XMEN_MSBIT);
    ln = (int)(n->len & XMEN_OBITS);
    if (ln == 1 && n->w[0] == 1)
        ln = 0;
    if ((ld + ln > xmen_mip->nib) && (xmen_mip->check || ld + ln > 2 * xmen_mip->nib))
        EasyXMen_berror(XMEN_ERR_FLASH_OVERFLOW);
    if (xmen_mip->ERNUM)
    {
        XMEN_OUT
        return;
    }
    copy(n, x);
    if (n->len == 0)
    {
        XMEN_OUT
        return;
    }
    s ^= (d->len & XMEN_MSBIT);
    if (ld == 0)
    {
        if (x->len != 0)
            x->len |= s;
        XMEN_OUT
        return;
    }
    for (i = 0; i < ld; i++)
        x->w[ln + i] = d->w[i];
    x->len = (s | (ln + ((unsigned int)ld << XMEN_BTS)));
    XMEN_OUT
}

/**
 * @description Changes the radix (base) of a multi-precision integer.
 *        This function converts the multi-precision integer `x` from an old base (`oldbase`) to the current base used by the system,
 *        storing the result in `y`. If the bases are the same, it simply copies `x` to `y`.
 *        The function handles both standard multi-precision integers and flash numbers, adjusting the representation as necessary.
 * @param[in] x         : Pointer to the multi-precision integer to change the radix of.
 * @param[in] oldbase   : The original base of the multi-precision integer `x`.
 * @param[out] y        : Pointer to the multi-precision integer where the converted value will be stored.
 * @return void
 */
static void cbase(XMen_BigTPtr x, XMen_SmallType oldbase, XMen_BigTPtr y)
{ /*  change radix of x from oldbase to base  */
    int i, s;
    XMen_SmallType n;
    BOOL done;

    if (xmen_mip->ERNUM)
        return;
    if (xmen_mip->base == oldbase)
    {
        copy(x, y);
        return;
    }

    XMEN_IN(13)

    s = exsign(x);
#ifdef XMEN_FLASH
    numer(x, xmen_mip->w1);
    denom(x, xmen_mip->w2);
    done = FALSE;
#else
    copy(x, xmen_mip->w1);
    done = TRUE;
#endif
    insign(PLUS, xmen_mip->w1);

    for (;;)
    {
        zero(xmen_mip->w6);
        convert(1, xmen_mip->w0);

        for (i = 0; i < (int)xmen_mip->w1->len; i++)
        { /* this is a bit slow - but not time critical */

            xmen_pmul(xmen_mip->w0, xmen_mip->w1->w[i], xmen_mip->w5);

            XMen_Add(xmen_mip->w6, xmen_mip->w5, xmen_mip->w6);
            if (oldbase == 0)
            { /* bit of a frig! */
                n = xmen_shiftbits(1, MIRACL / 2);
                xmen_pmul(xmen_mip->w0, n, xmen_mip->w0);
                xmen_pmul(xmen_mip->w0, n, xmen_mip->w0);
            }
            else
                xmen_pmul(xmen_mip->w0, oldbase, xmen_mip->w0);
        }
        if (xmen_mip->ERNUM || done)
            break;
#ifdef XMEN_FLASH
        copy(xmen_mip->w2, xmen_mip->w1);
        copy(xmen_mip->w6, xmen_mip->w7);
        done = TRUE;
#endif
    }

#ifdef XMEN_FLASH
    fpack(xmen_mip->w7, xmen_mip->w6, y);
#else
    copy(xmen_mip->w6, y);
#endif

    insign(s, y);
    XMEN_OUT
}

/**
 * @description Extracts a packed digit from a multi-precision integer.
 *        This function retrieves the `i`th digit from the multi-precision integer `x`, where digits are packed into words according to the system's packing scheme.
 *        The digit is returned as an integer value. If the packing factor (`xmen_mip->pack`) is 1, it directly returns the word value as the digit.
 *        Otherwise, it divides the word by the appropriate base power to isolate the desired digit.
 * @param[in] x   : Pointer to the multi-precision integer from which to extract the digit.
 * @param[in] i   : The position of the digit to extract (1-based index).
 * @return int    : The extracted digit as an integer.
 * @retval -      : Returns the extracted digit as an integer.
 */
static int getdig(XMen_BigTPtr x, int i)
{ /* extract a packed digit */
    int k;
    XMen_SmallType n;

    i--;
    n = x->w[i / xmen_mip->pack];

    if (xmen_mip->pack == 1)
        return (int)n;
    k = i % xmen_mip->pack;
    for (i = 1; i <= k; i++)
        n = XMEN_DIV(n, xmen_mip->apbase);
    return (int)XMEN_REMAIN(n, xmen_mip->apbase);
}

/**
 * @description Inserts a digit into a packed word of a multi-precision integer.
 *        This function sets the `i`th digit of the multi-precision integer `x` to the given integer `n`. Digits are packed into words according to the system's packing scheme.
 *        It ensures that the new digit is correctly placed within the appropriate word and adjusts the length of `x` if necessary.
 *        It also checks for overflow conditions and reports errors accordingly.
 * @param[in] n   : The integer value of the digit to insert.
 * @param[inout] x: Pointer to the multi-precision integer where the digit will be inserted.
 * @param[in] i   : The position of the digit to insert (1-based index).
 * @return void
 */
static void putdig(int n, XMen_BigTPtr x, int i)
{ /* insert a digit into a packed word */
    int j, k, lx;
    XMen_SmallType m, p;
    unsigned int s;

    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(26)

    s = (x->len & (XMEN_MSBIT));
    lx = (int)(x->len & (XMEN_OBITS));
    m = getdig(x, i);
    p = n;
    i--;
    j = i / xmen_mip->pack;
    k = i % xmen_mip->pack;
    for (i = 1; i <= k; i++)
    {
        m *= xmen_mip->apbase;
        p *= xmen_mip->apbase;
    }
    if (j >= xmen_mip->nib && (xmen_mip->check || j >= 2 * xmen_mip->nib))
    {
        EasyXMen_berror(XMEN_ERR_OVERFLOW);
        XMEN_OUT
        return;
    }

    x->w[j] = (x->w[j] - m) + p;
    if (j >= lx)
        x->len = ((j + 1) | s);
    xmen_lzero(x);
    XMEN_OUT
}

/**
 * @description Parses and inputs a big number from a string.
 *        This function reads a multi-precision integer or flash number from the provided string `string` and stores it in the flash number `x`.
 *        It handles various input formats, including ASCII and binary representations, and supports different numerical bases.
 *        The function returns the length of the parsed input in digits or bytes.
 * @param[inout] x      : Pointer to the flash number where the parsed value will be stored.
 * @param[in] string    : Pointer to the null-terminated string containing the number to parse.
 * @return int          : The length of the parsed input in digits or bytes.
 * @retval 0            : If an error occurs during parsing.
 */
static int instr(flash x, char *string)
{ /*  input a big number       *
   *  returns length in digits */
    int i, ipt, n, s, e;
    int ch;
#ifdef XMEN_FLASH
    BOOL frac;
#endif

    if (xmen_mip->ERNUM)
        return 0;

    XMEN_IN(76)

    if (xmen_mip->apbase == 0 || xmen_mip->apbase > 256)
    {
        EasyXMen_berror(XMEN_ERR_BASE_TOO_BIG);
        XMEN_OUT
        return 0;
    }

    if (!xmen_mip->active)
    {
        EasyXMen_berror(XMEN_ERR_NO_MIRSYS);
        XMEN_OUT
        return 0;
    }

    zero(x);
    if (xmen_mip->fin)
        string = xmen_mip->IOBUFF;
    if (xmen_mip->INPLEN == 0)
    { /* inputting ASCII bytes */
        for (;;)
        { /*  get input length  */
            ch = (unsigned char)string[xmen_mip->INPLEN];
            if (ch == '\0')
                break;
            if (xmen_mip->apbase <= 60 && ch == '\n')
                break;
            xmen_mip->INPLEN++;
            if (string == xmen_mip->IOBUFF && xmen_mip->INPLEN >= xmen_mip->IOBSIZ)
            {
                EasyXMen_berror(XMEN_ERR_IO_OVERFLOW);
                XMEN_OUT
                return 0;
            }
        }
    }
    else
    { /* inputting BINARY bytes */
        if (string == xmen_mip->IOBUFF && xmen_mip->INPLEN >= xmen_mip->IOBSIZ)
        {
            EasyXMen_berror(XMEN_ERR_IO_OVERFLOW);
            XMEN_OUT
            return 0;
        }
    }
    n = 0;
    s = PLUS;
    e = 0;
#ifdef XMEN_FLASH
    frac = FALSE;
#endif
    if (xmen_mip->INPLEN > 0 && xmen_mip->apbase <= 60)
    { /* skip leading blanks and check sign */
#ifdef XMEN_FLASH
        if (string[xmen_mip->INPLEN - 1] == '/')
            xmen_mip->INPLEN--;
#endif
        while (string[e] == ' ')
            e++;
        if (string[e] == '-')
        { /* check sign */
            s = MINUS;
            e++;
        }
        if (string[e] == '+')
            e++;
    }
    for (i = xmen_mip->INPLEN - 1; i >= e; i--)
    {
        ch = (unsigned char)string[i];
        if (xmen_mip->apbase <= 60 || xmen_mip->apbase == 64)
        { /* check for slash or dot and convert character to number */
#ifdef XMEN_FLASH
            if (xmen_mip->apbase <= 60 && !frac)
            {
                if (ch == '/')
                {
                    frac = TRUE;
                    copy(x, xmen_mip->w0);
                    zero(x);
                    n = 0;
                    continue;
                }
                if (ch == '.')
                {
                    frac = TRUE;
                    zero(xmen_mip->w0);
                    putdig(1, xmen_mip->w0, n + 1);
                    continue;
                }
            }
#endif
            ch += 80;
            if (xmen_mip->apbase == 64)
            { /* base64 */
                if (ch <= 112)
                    continue; /* ignore white space */
                if (ch > 144 && ch < 171)
                    ch -= 145;
                if (ch > 176 && ch < 203)
                    ch -= 151;
                if (ch > 127 && ch < 138)
                    ch -= 76;
                if (ch == 123)
                    ch = 62;
                if (ch == 127)
                    ch = 63;
                if (ch == 141)
                    continue; /* ignore pads '=' */
            }
            else
            {
                if (ch > 127 && ch < 138)
                    ch -= 128;
                if (ch > 144 && ch < 171)
                    ch -= 135;
                if (xmen_mip->apbase <= 16)
                {
                    if (ch > 176 && ch < 203)
                        ch -= 167;
                }
                else
                {
                    if (ch > 176 && ch < 203)
                        ch -= 141;
                }
            }
        }

        if ((XMen_SmallType)ch >= xmen_mip->apbase)
        {
            EasyXMen_berror(XMEN_ERR_BAD_FORMAT);
            XMEN_OUT
            return 0;
        }
        n++;
        putdig(ch, x, n);
    }
    ipt = xmen_mip->INPLEN;
    xmen_mip->INPLEN = 0;
    insign(s, x);
    xmen_lzero(x);
#ifdef XMEN_FLASH
    xmen_lzero(xmen_mip->w0);
    if (frac)
        fpack(x, xmen_mip->w0, x);
#endif
    XMEN_OUT
    return ipt;
}

/**
 * @description Strips leading zeros from a multi-precision integer.
 *        This function removes any leading zero words from the multi-precision integer `x` and adjusts its length accordingly.
 *        It preserves the sign bit of the number, ensuring that the sign is not altered during this process.
 * @param[inout] x : Pointer to the multi-precision integer from which to remove leading zeros.
 * @return void
 */
static void xmen_lzero(XMen_BigTPtr x)
{ /*  strip leading zeros from big number  */
    unsigned int s;
    int m;
    s = (x->len & (XMEN_MSBIT));
    m = (int)(x->len & (XMEN_OBITS));
    while (m > 0 && x->w[m - 1] == 0)
        m--;
    x->len = m;
    if (m > 0)
        x->len |= s;
}

/**
 * @description Performs division of a multi-precision integer by a small integer.
 *        This function divides the multi-precision integer `x` by the small integer `sn`, storing the quotient in `z` and returning the remainder.
 *        Depending on the configuration, it can also handle full-width operations and floating-point rounding for more precise results.
 * @param[in] x   : Pointer to the multi-precision integer dividend.
 * @param[in] sn  : The small integer divisor.
 #ifdef XMEN_FP_ROUNDING
 * @param[in] isn : The inverse of the small integer divisor, used for floating-point rounding.
 #endif
 * @param[out] z  : Pointer to the multi-precision integer where the quotient will be stored.
 * @return XMen_SmallType : The remainder of the division.
 */
#ifdef XMEN_FP_ROUNDING
static XMen_SmallType xmen_sdiv(XMen_BigTPtr x, XMen_SmallType sn, XMen_LargeType isn, XMen_BigTPtr z)
#else
static XMen_SmallType xmen_sdiv(XMen_BigTPtr x, XMen_SmallType sn, XMen_BigTPtr z)
#endif
{
    int i, xl;
    XMen_SmallType sr, *xg, *zg;
    union doubleword dble;
    XMen_LargeType dbled;
    XMen_LargeType ldres;

    sr = 0;
    xl = (int)(x->len & XMEN_OBITS);
    if (x != z)
        zero(z);
#ifndef XMEN_SIMPLE_BASE
    if (xmen_mip->base == 0)
    {
#endif
#ifndef XMEN_NOFULLWIDTH
        xg = x->w;
        zg = z->w;
        /* inline - substitutes for loop below */

        for (i = xl - 1; i >= 0; i--)
        {
            dble.h[BYTE_BOT] = x->w[i];
            dble.h[BYTE_TOP] = sr;
            z->w[i] = (XMen_SmallType)(dble.d / sn);
            sr = (XMen_SmallType)(dble.d - (XMen_LargeType)z->w[i] * sn);
        }
#endif
#ifndef XMEN_SIMPLE_BASE
    }
    else
        for (i = xl - 1; i >= 0; i--)
        { /* divide each digit of x by n */
            dbled = (XMen_LargeType)sr * xmen_mip->base + x->w[i];
#ifdef XMEN_FP_ROUNDING
            z->w[i] = (XMen_SmallType)XMEN_LROUND(dbled * isn);
#else
            z->w[i] = (XMen_SmallType)XMEN_LROUND(dbled / sn);
#endif
            sr = (XMen_SmallType)(dbled - (XMen_LargeType)z->w[i] * sn);
        }
#endif
    z->len = x->len;
    xmen_lzero(z);
    return sr;
}

/**
 * @description Normalizes the divisor for efficient multi-precision division.
 *        This function normalizes the multi-precision integer `y` by scaling it such that its most significant digit is close to the base value.
 *        It optionally copies the value of `x` to `y` if they are not the same, and then calculates a normalization factor.
 *        If the normalization factor is not 1, it multiplies `y` by this factor to achieve the normalized form.
 * @param[in] x   : Pointer to the multi-precision integer from which to copy (if different from `y`).
 * @param[inout] y: Pointer to the multi-precision integer to normalize.
 * @return XMen_SmallType : The normalization factor used.
 */
static XMen_SmallType normalise(XMen_BigTPtr x, XMen_BigTPtr y)
{ /* normalise divisor */
    XMen_SmallType norm, r;
    int len;

    XMEN_IN(4)

    if (x != y)
        copy(x, y);
    len = (int)(y->len & XMEN_OBITS);
#ifndef XMEN_SIMPLE_BASE
    if (xmen_mip->base == 0)
    {
#endif
#ifndef XMEN_NOFULLWIDTH
        if ((r = y->w[len - 1] + 1) == 0)
            norm = 1;
        else
            norm = (XMen_SmallType)(((XMen_LargeType)1 << MIRACL) / r);
        if (norm != 1)
            xmen_pmul(y, norm, y);
#endif
#ifndef XMEN_SIMPLE_BASE
    }
    else
    {
        norm = XMEN_DIV(xmen_mip->base, (XMen_SmallType)(y->w[len - 1] + 1));
        if (norm != 1)
            xmen_pmul(y, norm, y);
    }
#endif
    XMEN_OUT
    return norm;
}

/**
 * @description Divides a multi-precision integer by an integer and returns the quotient and remainder.
 *        This function divides the multi-precision integer `x` by the integer `n`, storing the quotient in `z` and returning the remainder.
 *        It handles special cases, such as division by 1 and fast division by 2 using shifting, and supports floating-point rounding for more precise results.
 * @param[in] x   : Pointer to the multi-precision integer dividend.
 * @param[in] n   : The integer divisor.
 * @param[out] z  : Pointer to the multi-precision integer where the quotient will be stored.
 * @return int    : The remainder of the division.
 * @retval -      : Returns the remainder as an integer. If `x` is negative, it returns the negated remainder.
 */
static int subdiv(XMen_BigTPtr x, int n, XMen_BigTPtr z)
{ /*  subdivide a big number by an int   z=x/n  *
   *  returns int remainder                     */
    unsigned int sx;
#ifdef XMEN_FP_ROUNDING
    XMen_LargeType in;
#endif
    int r, i, msb;
    XMen_SmallType lsb;

    if (xmen_mip->ERNUM)
        return 0;

    XMEN_IN(10)
#ifdef XMEN_FLASH
    if (xmen_notint(x))
        EasyXMen_berror(XMEN_ERR_INT_OP);
#endif
    if (n == 0)
        EasyXMen_berror(XMEN_ERR_DIV_BY_ZERO);
    if (xmen_mip->ERNUM)
    {
        XMEN_OUT
        return 0;
    }

    if (x->len == 0)
    {
        zero(z);
        XMEN_OUT
        return 0;
    }
    if (n == 1) /* special case */
    {
        copy(x, z);
        XMEN_OUT
        return 0;
    }
    sx = (x->len & XMEN_MSBIT);
    if (n == 2 && xmen_mip->base == 0)
    { /* fast division by 2 using shifting */
#ifndef XMEN_NOFULLWIDTH

        /* I don't want this code upsetting the compiler ... */
        /* xmen_mip->base==0 can't happen with XMEN_NOFULLWIDTH  */

        copy(x, z);
        msb = (int)(z->len & XMEN_OBITS) - 1;
        r = (int)z->w[0] & 1;
        for (i = 0;; i++)
        {
            z->w[i] >>= 1;
            if (i == msb)
            {
                if (z->w[i] == 0)
                    xmen_lzero(z);
                break;
            }
            lsb = z->w[i + 1] & 1;
            z->w[i] |= (lsb << (MIRACL - 1));
        }

        XMEN_OUT
        if (sx == 0)
            return r;
        else
            return (-r);
#endif
    }

#ifdef XMEN_FP_ROUNDING
    in = xmen_invert(n);
#endif
    if (n < 0)
    {
        n = (-n);
#ifdef XMEN_FP_ROUNDING
        r = (int)xmen_sdiv(x, (XMen_SmallType)n, in, z);
#else
        r = (int)xmen_sdiv(x, (XMen_SmallType)n, z);
#endif
        if (z->len != 0)
            z->len ^= XMEN_MSBIT;
    }
#ifdef XMEN_FP_ROUNDING
    else
        r = (int)xmen_sdiv(x, (XMen_SmallType)n, in, z);
#else
    else
        r = (int)xmen_sdiv(x, (XMen_SmallType)n, z);
#endif
    XMEN_OUT
    if (sx == 0)
        return r;
    else
        return (-r);
}

/**
 * @description Returns the integer remainder of a multi-precision integer divided by an integer.
 *        This function calculates the remainder when the multi-precision integer `x` is divided by the integer `n`.
 *        It includes optimizations for specific divisors like 2 and 8, using bitwise operations for faster computation.
 * @param[in] x   : Pointer to the multi-precision integer dividend.
 * @param[in] n   : The integer divisor.
 * @return int    : The integer remainder of the division.
 * @retval -      : Returns the remainder as an integer. If `x` is negative, it returns the negated remainder.
 */
static int remain(XMen_BigTPtr x, int n)
{ /* return integer remainder when x divided by n */
    int r;
    unsigned int sx;

    if (xmen_mip->ERNUM)
        return FALSE;

    XMEN_IN(88);

    sx = (x->len & XMEN_MSBIT);

    if (n == 2 && XMEN_REMAIN(xmen_mip->base, 2) == 0)
    { /* fast odd/even check if base is even */
        XMEN_OUT
        if ((int)XMEN_REMAIN(x->w[0], 2) == 0)
            return 0;
        else
        {
            if (sx == 0)
                return 1;
            else
                return (-1);
        }
    }
    if (n == 8 && XMEN_REMAIN(xmen_mip->base, 8) == 0)
    { /* fast check */
        XMEN_OUT
        r = (int)XMEN_REMAIN(x->w[0], 8);
        if (sx != 0)
            r = -r;
        return r;
    }

    copy(x, xmen_mip->w0);
    r = subdiv(xmen_mip->w0, n, xmen_mip->w0);
    XMEN_OUT
    return r;
}

/**
 * @description Performs a multiplication followed by division with a remainder.
 *        This function multiplies two small integers `a` and `c`, then divides the result by `m`, returning the quotient and storing the remainder in `rp`.
 * @param[in] a   : The first small integer to multiply.
 * @param[in] c   : The second small integer to multiply.
 * @param[in] m   : The divisor used after multiplication.
 * @param[out] rp : Pointer to store the remainder of the division.
 * @return XMen_SmallType : The quotient of the division.
 */
static XMen_SmallType muldvm(XMen_SmallType a, XMen_SmallType c, XMen_SmallType m, XMen_SmallType *rp)
{
    XMen_SmallType q;
    union doubleword dble;
    dble.h[BYTE_BOT] = c;
    dble.h[BYTE_TOP] = a;

    q = (XMen_SmallType)(dble.d / m);
    *rp = (XMen_SmallType)(dble.d - (XMen_LargeType)q * m);
    return q;
}

/**
 * @description Performs a fast division of two large integers, expecting a small quotient.
 *        This function divides the large integer `u` by the large integer `v` and returns the quotient.
 *        It optimizes for small quotients by subtracting `v` from `u` up to eight times before falling back to a more complex division algorithm.
 * @param[in] u   : The large integer dividend.
 * @param[in] v   : The large integer divisor.
 * @return XMen_SmallType : The quotient of the division.
 */
static XMen_SmallType qdiv(XMen_LargeType u, XMen_LargeType v)
{ /* fast division - small quotient expected.  */
    XMen_LargeType lq, x = u;
    x -= v;
    if (x < v)
        return 1;
    x -= v;
    if (x < v)
        return 2;
    x -= v;
    if (x < v)
        return 3;
    x -= v;
    if (x < v)
        return 4;
    x -= v;
    if (x < v)
        return 5;
    x -= v;
    if (x < v)
        return 6;
    x -= v;
    if (x < v)
        return 7;
    x -= v;
    if (x < v)
        return 8;

    /* do it the hard way! */

    lq = 8 + XMEN_DIV(x, v);
    if (lq >= MAXBASE)
        return 0;
    return (XMen_SmallType)lq;
}

/**
 * @description Negates a multi-precision or flash number.
 *        This function negates the multi-precision or flash number `x`, storing the result in `y`.
 *        It simply copies `x` to `y` and toggles the sign bit if `y` is non-zero.
 * @param[in] x   : Pointer to the multi-precision or flash number to negate.
 * @param[out] y  : Pointer to the multi-precision or flash number where the negated value will be stored.
 * @return void
 */
static void negify(flash x, flash y)
{ /* negate a big/flash variable: y=-x */
    copy(x, y);
    if (y->len != 0)
        y->len ^= XMEN_MSBIT;
}

/**
 * @description Computes the modular multiplicative inverse of a multi-precision integer.
 *        This function calculates the modular multiplicative inverse of `x` modulo `y`, storing the result in `z`.
 *        It returns the greatest common divisor (GCD) of `x` and `y`.
 * @param[in] x   : Pointer to the multi-precision integer for which to find the modular inverse.
 * @param[in] y   : Pointer to the modulus.
 * @param[out] z  : Pointer to the multi-precision integer where the result will be stored.
 * @return int    : The greatest common divisor (GCD) of `x` and `y`.
 * @retval -      : Returns the GCD; if the GCD is 1, `z` contains the modular inverse.
 */
static int invmodp(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr z)
{

    int gcd;

    XMEN_IN(213);
    gcd = XMen_Xgcd(x, y, z, z, z);
    XMEN_OUT
    return gcd;
}

/**
 * @description Shifts a multi-precision integer by multiplying it with `xmen_base^n`.
 *        This function shifts the multi-precision integer `x` by `n` places, effectively multiplying it by `xmen_base^n` and stores the result in `w`.
 *        Positive `n` values shift left (multiplication), while negative `n` values shift right (division).
 * @param[in] x   : Pointer to the multi-precision integer to shift.
 * @param[in] n   : The number of places to shift. Positive for left shift, negative for right shift.
 * @param[out] w  : Pointer to the multi-precision integer where the shifted value will be stored.
 * @return void
 */
static void xmen_shift(XMen_BigTPtr x, int n, XMen_BigTPtr w)
{ /* set w=x.(xmen_base^n) by shifting */
    unsigned int s;
    int i, bl;
    XMen_SmallType *gw = w->w;

    if (xmen_mip->ERNUM)
        return;
    copy(x, w);
    if (w->len == 0 || n == 0)
        return;
    XMEN_IN(33)

    if (xmen_notint(w))
        EasyXMen_berror(XMEN_ERR_INT_OP);
    s = (w->len & (XMEN_MSBIT));
    bl = (int)(w->len & (XMEN_OBITS)) + n;
    if (bl <= 0)
    {
        zero(w);
        XMEN_OUT
        return;
    }
    if (bl > xmen_mip->nib && xmen_mip->check)
        EasyXMen_berror(XMEN_ERR_OVERFLOW);
    if (xmen_mip->ERNUM)
    {
        XMEN_OUT
        return;
    }
    if (n > 0)
    {
        for (i = bl - 1; i >= n; i--)
            gw[i] = gw[i - n];
        for (i = 0; i < n; i++)
            gw[i] = 0;
    }
    else
    {
        n = (-n);
        for (i = 0; i < bl; i++)
            gw[i] = gw[i + n];
        for (i = 0; i < n; i++)
            gw[bl + i] = 0;
    }
    w->len = (bl | s);
    XMEN_OUT
}

/**
 * @description Converts a multi-precision integer to n-residue format.
 *        This function converts the multi-precision integer `x` into n-residue format, storing the result in `y`.
 *        It ensures that the modulus is set and handles the conversion process, including normalization for Montgomery multiplication if enabled.
 * @param[in] x   : Pointer to the multi-precision integer to convert.
 * @param[out] y  : Pointer to the multi-precision integer where the converted value will be stored.
 * @return void
 */
static void nres(XMen_BigTPtr x, XMen_BigTPtr y)
{ /* convert x to n-residue format */

    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(81)

    if (EasyXMen_sm2_size(xmen_mip->modulus) == 0)
    {
        EasyXMen_berror(XMEN_ERR_NO_MODULUS);
        XMEN_OUT
        return;
    }
    copy(x, y);
    XMen_Divide(y, xmen_mip->modulus, xmen_mip->modulus);
    if (EasyXMen_sm2_size(y) < 0)
        XMen_Add(y, xmen_mip->modulus, y);
    if (!xmen_mip->MONTY)
    {
        XMEN_OUT
        return;
    }
    xmen_mip->check = OFF;

    xmen_shift(y, (int)xmen_mip->modulus->len, xmen_mip->w0);
    XMen_Divide(xmen_mip->w0, xmen_mip->modulus, xmen_mip->modulus);
    xmen_mip->check = ON;
    copy(xmen_mip->w0, y);

    XMEN_OUT
}

/* prepare Montgomery modulus */
/**
 * @description Prepares the system for Montgomery multiplication with a given modulus.
 *        This function sets up the necessary parameters for performing Montgomery arithmetic operations using the provided modulus `n`.
 *        It initializes various internal states and checks, ensuring that the modulus is suitable for use in cryptographic operations.
 * @param[in] n   : Pointer to the multi-precision integer modulus to prepare for Montgomery arithmetic.
 * @return XMen_SmallType : Returns the value of `xmen_mip->ndash`, which is used in Montgomery reduction.
 * @retval 0      : If an error occurs during setup or if Montgomery arithmetic is disabled.
 */
XMen_SmallType prepare_monty(XMen_BigTPtr n)
{
    if (xmen_mip->ERNUM)
        return (XMen_SmallType)0;
    /* Is it set-up already? */
    if (EasyXMen_sm2_size(xmen_mip->modulus) != 0)
        if (XMen_Compare(n, xmen_mip->modulus) == 0)
            return xmen_mip->ndash;

    XMEN_IN(80)

    if (EasyXMen_sm2_size(n) <= 2)
    {
        EasyXMen_berror(XMEN_ERR_BAD_MODULUS);
        XMEN_OUT
        return (XMen_SmallType)0;
    }

    zero(xmen_mip->w6);
    zero(xmen_mip->w15);

    /* set a small negative QNR (on the assumption that n is prime!) */
    /* These defaults can be over-ridden                             */

    /* Did you know that for p=2 mod 3, -3 is a QNR? */

    xmen_mip->pmod8 = remain(n, 8);

    switch (xmen_mip->pmod8)
    {
    case 0:
    case 1:
    case 2:
    case 4:
    case 6:
        xmen_mip->qnr = 0; /* none defined */
        break;
    case 3:
        xmen_mip->qnr = -1;
        break;
    case 5:
        xmen_mip->qnr = -2;
        break;
    case 7:
        xmen_mip->qnr = -1;
        break;
    }
    xmen_mip->pmod9 = remain(n, 9);

    xmen_mip->NO_CARRY = FALSE;
    if (n->w[n->len - 1] >> (MIRACL - 4) < 5)
        xmen_mip->NO_CARRY = TRUE;

#ifdef XMEN_DISABLE_MONTGOMERY
    xmen_mip->MONTY = OFF;
#else
    xmen_mip->MONTY = ON;
#endif

#ifdef XMEN_COMBA
    xmen_mip->ACTIVE = FALSE;

    if (XMEN_COMBA == n->len && xmen_mip->base == xmen_mip->base2)
    {
        xmen_mip->ACTIVE = TRUE;
#ifdef XMEN_SPECIAL
        xmen_mip->MONTY = OFF; /* "special" modulus reduction */

#endif /* implemented in mrcomba.c    */
    }

#endif
    convert(1, xmen_mip->one);
    if (!xmen_mip->MONTY)
    { /* Montgomery arithmetic is turned off */
        copy(n, xmen_mip->modulus);
        xmen_mip->ndash = 0;
        XMEN_OUT
        return (XMen_SmallType)0;
    }
    xmen_mip->w6->len = 2;
    xmen_mip->w6->w[0] = 0;
    xmen_mip->w6->w[1] = 1; /* w6 = base */
    xmen_mip->w15->len = 1;
    xmen_mip->w15->w[0] = n->w[0]; /* w15 = n mod base */
    if (invmodp(xmen_mip->w15, xmen_mip->w6, xmen_mip->w14) != 1)
    { /* problems */
        EasyXMen_berror(XMEN_ERR_BAD_MODULUS);
        XMEN_OUT
        return (XMen_SmallType)0;
    }

    xmen_mip->ndash = xmen_mip->base - xmen_mip->w14->w[0]; /* = N' mod b */
    copy(n, xmen_mip->modulus);
    xmen_mip->check = OFF;
    xmen_shift(xmen_mip->modulus, (int)xmen_mip->modulus->len, xmen_mip->pR);
    xmen_mip->check = ON;
    nres(xmen_mip->one, xmen_mip->one);
    XMEN_OUT

    return xmen_mip->ndash;
}

#ifdef XMEN_FP_ROUNDING
/**
 * @description Performs a combined multiplication, addition, division, and remainder operation.
 *        This function calculates the quotient and remainder of ((a * b) + c) divided by m.
 *        It uses an inverse multiplier `im` for efficient division.
 * @param[in] a   : The first small integer to multiply.
 * @param[in] b   : The second small integer to multiply.
 * @param[in] c   : The integer to add to the product of `a` and `b`.
 * @param[in] m   : The divisor used after multiplication and addition.
 * @param[in] im  : The inverse multiplier used for efficient division.
 * @param[out] rp : Pointer to store the remainder of the division.
 * @return XMen_SmallType : The quotient of the division.
 */
static XMen_SmallType imuldiv(XMen_SmallType a, XMen_SmallType b, XMen_SmallType c, XMen_SmallType m, XMen_LargeType im, XMen_SmallType *rp)
{
    XMen_SmallType q;
    XMen_LargeType ldres, p = (XMen_LargeType)a * b + c;
    q = (XMen_SmallType)XMEN_LROUND(p * im);
    *rp = (XMen_SmallType)(p - (XMen_LargeType)q * m);
    return q;
}

#endif
/**
 * @description Performs Montgomery's REDC (Reduction) function to reduce a number modulo another.
 *        This function implements the REDC algorithm, which is used for efficient modular reduction in Montgomery arithmetic.
 *        It can also convert n-residues back to normal form. The function ensures that the result is within the allowed modulus range.
 * @param[in] x   : Pointer to the multi-precision integer to be reduced.
 * @param[out] y  : Pointer to the multi-precision integer where the reduced value will be stored.
 * @return void
 */
static void redc(XMen_BigTPtr x, XMen_BigTPtr y)
{ /* Montgomery's REDC function p. 520 */
    /* also used to convert n-residues back to normal form */
    XMen_SmallType carry, delay_carry, m, ndash, *w0g, *mg;

#ifdef XMEN_ITANIUM
    XMen_SmallType tm;
#endif
#ifdef XMEN_WIN64
    XMen_SmallType tm, tr;
#endif
    int i, j, rn, rn2;
    XMen_BigTPtr w0, modulus;
    union doubleword dble;
    XMen_LargeType dbled, ldres;

    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(82)

    w0 = xmen_mip->w0; /* get these into local variables (for inline assembly) */
    modulus = xmen_mip->modulus;
    ndash = xmen_mip->ndash;

    copy(x, w0);
    if (!xmen_mip->MONTY)
    {
        /*#ifdef XMEN_CELL
                mod256(w0,w0);
        #else */
        XMen_Divide(w0, modulus, modulus);
        /* #endif */
        copy(w0, y);
        XMEN_OUT
        return;
    }
    delay_carry = 0;
    rn = (int)modulus->len;
    rn2 = rn + rn;
#ifndef XMEN_SIMPLE_BASE
    if (xmen_mip->base == 0)
    {
#endif
#ifndef XMEN_NOFULLWIDTH
        mg = modulus->w;
        w0g = w0->w;
        for (i = 0; i < rn; i++)
        {
            /*        muldvd(w0->w[i],ndash,0,&m);    Note that after this time   */
            m = ndash * w0->w[i];
            carry = 0; /* around the loop, w0[i]=0    */

            for (j = 0; j < rn; j++)
            {
                dble.d = (XMen_LargeType)m * modulus->w[j] + carry + w0->w[i + j];
                w0->w[i + j] = dble.h[BYTE_BOT];
                carry = dble.h[BYTE_TOP];
            }
            w0->w[rn + i] += delay_carry;
            if (w0->w[rn + i] < delay_carry)
                delay_carry = 1;
            else
                delay_carry = 0;
            w0->w[rn + i] += carry;
            if (w0->w[rn + i] < carry)
                delay_carry = 1;
        }
#endif

#ifndef XMEN_SIMPLE_BASE
    }
    else
        for (i = 0; i < rn; i++)
        {
#ifdef XMEN_FP_ROUNDING
            imuldiv(w0->w[i], ndash, 0, xmen_mip->base, xmen_mip->inverse_base, &m);
#else
            muldiv(w0->w[i], ndash, 0, xmen_mip->base, &m);
#endif
            carry = 0;
            for (j = 0; j < rn; j++)
            {
                dbled = (XMen_LargeType)m * modulus->w[j] + carry + w0->w[i + j];
#ifdef XMEN_FP_ROUNDING
                carry = (XMen_SmallType)XMEN_LROUND(dbled * xmen_mip->inverse_base);
#else
                if (xmen_mip->base == xmen_mip->base2)
                    carry = (XMen_SmallType)(dbled >> xmen_mip->lg2b);
                else
                    carry = (XMen_SmallType)XMEN_LROUND(dbled / xmen_mip->base);
#endif
                w0->w[i + j] = (XMen_SmallType)(dbled - (XMen_LargeType)carry * xmen_mip->base);
            }
            w0->w[rn + i] += (delay_carry + carry);
            delay_carry = 0;
            if (w0->w[rn + i] >= xmen_mip->base)
            {
                w0->w[rn + i] -= xmen_mip->base;
                delay_carry = 1;
            }
        }
#endif
    w0->w[rn2] = delay_carry;
    w0->len = rn2 + 1;
    xmen_shift(w0, (-rn), w0);
    xmen_lzero(w0);

    if (XMen_Compare(w0, modulus) >= 0)
        xmen_psub(w0, modulus, w0);
    copy(w0, y);
    XMEN_OUT
}

/**
 * @description Performs modular multiplication using n-residues.
 *        This function calculates the product of two multi-precision integers `x` and `y`, modulo `n`, and stores the result in `w`.
 *        It utilizes n-residue arithmetic for efficient computation, which is particularly useful in cryptographic operations.
 * @param[in] x   : Pointer to the first multi-precision integer operand.
 * @param[in] y   : Pointer to the second multi-precision integer operand.
 * @param[out] w  : Pointer to the multi-precision integer where the result (x * y mod n) will be stored.
 * @return void
 */
static void nres_modmult(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr w)
{ /* Modular multiplication using n-residues w=x*y mod n */

    if ((x == NULL || x->len == 0) && x == w)
        return;
    if ((y == NULL || y->len == 0) && y == w)
        return;
    if (y == NULL || x == NULL || x->len == 0 || y->len == 0)
    {
        zero(w);
        return;
    }
#ifdef XMEN_COUNT_OPS
    fpc++;
#endif
#ifdef XMEN_COMBA
    if (xmen_mip->ACTIVE)
    {
        if (x == y)
            comba_square(x, xmen_mip->w0);
        else
            comba_mult(x, y, xmen_mip->w0);
        comba_redc(xmen_mip->w0, w);
    }
    else
    {
#endif
        if (xmen_mip->ERNUM)
            return;

        XMEN_IN(83)

        xmen_mip->check = OFF;
        XMen_Multiply(x, y, xmen_mip->w0);
        redc(xmen_mip->w0, w);
        xmen_mip->check = ON;
        XMEN_OUT
#ifdef XMEN_COMBA
    }
#endif
}

/**
 * @description Computes the right-hand side of the elliptic curve equation for a given x-coordinate.
 *        This function calculates $ y^2 = x^3 + Ax + B $ for an elliptic curve, where `x` is the input coordinate and `y` is used as temporary storage.
 *        The result is stored in `y`.
 * @param[in] x   : Pointer to the multi-precision integer representing the x-coordinate.
 * @param[out] y  : Pointer to the multi-precision integer where the result (right-hand side of the elliptic curve equation) will be stored.
 * @return void
 */
static void epoint_getrhs(XMen_BigTPtr x, XMen_BigTPtr y)
{ /* x and y must be different */

    /* find x^3+Ax+B */

    nres_modmult(x, x, y);

    nres_modmult(y, x, y);
    if (xmen_abs(xmen_mip->Asize) == XMEN_TOOBIG)
        nres_modmult(x, xmen_mip->A, xmen_mip->w1);
    else
        nres_premult(x, xmen_mip->Asize, xmen_mip->w1);
    nres_modadd(y, xmen_mip->w1, y);
    if (xmen_abs(xmen_mip->Bsize) == XMEN_TOOBIG)
        nres_modadd(y, xmen_mip->B, y);
    else
    {
        convert(xmen_mip->Bsize, xmen_mip->w1);
        nres(xmen_mip->w1, xmen_mip->w1);
        nres_modadd(y, xmen_mip->w1, y);
    }
}

/**
 * @description Subtracts an integer from a multi-precision integer.
 *        This function subtracts an integer `n` from the multi-precision integer `x` and stores the result in `z`.
 * @param[in] x   : Pointer to the multi-precision integer from which to subtract.
 * @param[in] n   : The integer value to subtract from `x`.
 * @param[out] z  : Pointer to the multi-precision integer where the result (x - n) will be stored.
 * @return void
 */
static void decr(XMen_BigTPtr x, int n, XMen_BigTPtr z)
{ /* subtract int from big number: z=x-n */

    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(8)

    convert(n, xmen_mip->w0);
    xmen_select(x, MINUS, xmen_mip->w0, z);

    XMEN_OUT
}

/**
 * @description Calculates the number of bits in a multi-precision integer.
 *        This function returns the bit length of the given multi-precision integer `x`.
 *        It handles both binary and non-binary bases for internal representation.
 * @param[in] x   : Pointer to the multi-precision integer whose bit length is to be calculated.
 * @return int    : The number of bits in the multi-precision integer `x`.
 * @retval 0      : If an error occurs or if `x` has a size of zero.
 */
static int logb2(XMen_BigTPtr x)
{ /* returns number of bits in x */
    int xl, lg2;
    XMen_SmallType top;

    if (xmen_mip->ERNUM || EasyXMen_sm2_size(x) == 0)
        return 0;

    XMEN_IN(49)

#ifndef XMEN_ALWAYS_BINARY
    if (xmen_mip->base == xmen_mip->base2)
    {
#endif
        xl = (int)(x->len & XMEN_OBITS);
        lg2 = xmen_mip->lg2b * (xl - 1);
        top = x->w[xl - 1];
        while (top >= 1)
        {
            lg2++;
            top /= 2;
        }

#ifndef XMEN_ALWAYS_BINARY
    }
    else
    {
        copy(x, xmen_mip->w0);
        insign(PLUS, xmen_mip->w0);
        lg2 = 0;
        while (xmen_mip->w0->len > 1)
        {
            xmen_sdiv(xmen_mip->w0, xmen_mip->base2, xmen_mip->w0);

            lg2 += xmen_mip->lg2b;
        }

        while (xmen_mip->w0->w[0] >= 1)
        {
            lg2++;
            xmen_mip->w0->w[0] /= 2;
        }
    }
#endif
    XMEN_OUT
    return lg2;
}

/**
 * @description Tests the value of the n-th bit in a multi-precision integer.
 *        This function checks if the n-th bit (0-indexed) of the multi-precision integer `x` is set.
 * @param[in] x   : Pointer to the multi-precision integer whose bit is to be tested.
 * @param[in] n   : The index of the bit to test.
 * @return int    : Returns 1 if the n-th bit is set, otherwise returns 0.
 */
static int xmen_testbit(XMen_BigTPtr x, int n)
{ /* return value of n-th bit of big */

    if ((x->w[n / xmen_mip->lg2b] & ((XMen_SmallType)1 << (n % xmen_mip->lg2b))) > 0)
        return 1;

    return 0;
}

/**
 * @description Performs modular subtraction on two multi-precision integers.
 *        This function calculates $ w = x - y \mod n $ and stores the result in `w`.
 *        If `x < y`, it adjusts the result to ensure it remains within the modulus range.
 * @param[in] x   : Pointer to the first multi-precision integer operand.
 * @param[in] y   : Pointer to the second multi-precision integer operand.
 * @param[out] w  : Pointer to the multi-precision integer where the result will be stored.
 * @return void
 */
static void nres_modsub(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr w)
{ /* modular subtraction */

#ifdef XMEN_COMBA
    if (xmen_mip->ACTIVE)
    {
        comba_modsub(x, y, w);
        return;
    }
    else
    {
#endif
        if (xmen_mip->ERNUM)
            return;

        XMEN_IN(91)

        if (XMen_Compare(x, y) >= 0)
            xmen_psub(x, y, w);
        else
        {
            xmen_psub(y, x, w);
            xmen_psub(xmen_mip->modulus, w, w);
        }

        XMEN_OUT
#ifdef XMEN_COMBA
    }
#endif
}

/**
 * @description Sets a multi-precision integer to $ 2^n $.
 *        This function calculates $ x = 2^n $ and stores the result in `x`.
 * @param[in] n   : The exponent to which 2 is raised.
 * @param[out] x  : Pointer to the multi-precision integer where the result $ 2^n $ will be stored.
 * @return void
 */
static void expb2(int n, XMen_BigTPtr x)
{ /* sets x=2^n */
    int r, p;
#ifndef XMEN_ALWAYS_BINARY
    int i;
#endif
    if (xmen_mip->ERNUM)
        return;
    convert(1, x);
    if (n == 0)
        return;

    XMEN_IN(149)

    if (n < 0)
    {
        EasyXMen_berror(XMEN_ERR_NEG_POWER);
        XMEN_OUT
        return;
    }
    r = n / xmen_mip->lg2b;
    p = n % xmen_mip->lg2b;

#ifndef XMEN_ALWAYS_BINARY
    if (xmen_mip->base == xmen_mip->base2)
    {
#endif
        xmen_shift(x, r, x);
        x->w[x->len - 1] = xmen_shiftbits(x->w[x->len - 1], p);
#ifndef XMEN_ALWAYS_BINARY
    }
    else
    {
        for (i = 1; i <= r; i++)
            xmen_pmul(x, xmen_mip->base2, x);
        xmen_pmul(x, xmen_shiftbits((XMen_SmallType)1, p), x);
    }
#endif
    XMEN_OUT
}

/**
 * @description Computes the Lucas sequence values for a given multi-precision integer.
 *        This function calculates the Lucas sequence $ V_r $ and $ V_{r-1} $ for a given exponent `r`.
 *        The results are stored in `v` and `vp` respectively.
 * @param[in] p   : Pointer to the multi-precision integer representing the initial value of the sequence.
 * @param[in] r   : Pointer to the multi-precision integer representing the exponent.
 * @param[out] vp : Pointer to the multi-precision integer where $ V_{r-1} $ will be stored.
 * @param[out] v  : Pointer to the multi-precision integer where $ V_r $ will be stored.
 * @return void
 */
static void nres_lucas(XMen_BigTPtr p, XMen_BigTPtr r, XMen_BigTPtr vp, XMen_BigTPtr v)
{
    int i, nb;
    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(107)

    if (EasyXMen_sm2_size(r) == 0)
    {
        zero(vp);
        convert(2, v);
        nres(v, v);
        XMEN_OUT
        return;
    }
    if (EasyXMen_sm2_size(r) == 1 || EasyXMen_sm2_size(r) == (-1))
    { /* note - sign of r doesn't matter */
        convert(2, vp);
        nres(vp, vp);
        copy(p, v);
        XMEN_OUT
        return;
    }

    copy(p, xmen_mip->w3);

    convert(2, xmen_mip->w4);
    nres(xmen_mip->w4, xmen_mip->w4); /* w4=2 */

    copy(xmen_mip->w4, xmen_mip->w8);
    copy(xmen_mip->w3, xmen_mip->w9);

    copy(r, xmen_mip->w1);
    insign(PLUS, xmen_mip->w1);
    decr(xmen_mip->w1, 1, xmen_mip->w1);

#ifndef XMEN_ALWAYS_BINARY
    if (xmen_mip->base == xmen_mip->base2)
    {
#endif
        nb = logb2(xmen_mip->w1);
        for (i = nb - 1; i >= 0; i--)
        {
            if (xmen_mip->user != NULL)
                (*xmen_mip->user)();

            if (xmen_testbit(xmen_mip->w1, i))
            {
                nres_modmult(xmen_mip->w8, xmen_mip->w9, xmen_mip->w8);
                nres_modsub(xmen_mip->w8, xmen_mip->w3, xmen_mip->w8);
                nres_modmult(xmen_mip->w9, xmen_mip->w9, xmen_mip->w9);
                nres_modsub(xmen_mip->w9, xmen_mip->w4, xmen_mip->w9);
            }
            else
            {
                nres_modmult(xmen_mip->w9, xmen_mip->w8, xmen_mip->w9);
                nres_modsub(xmen_mip->w9, xmen_mip->w3, xmen_mip->w9);
                nres_modmult(xmen_mip->w8, xmen_mip->w8, xmen_mip->w8);
                nres_modsub(xmen_mip->w8, xmen_mip->w4, xmen_mip->w8);
            }
        }

#ifndef XMEN_ALWAYS_BINARY
    }
    else
    {
        expb2(logb2(xmen_mip->w1) - 1, xmen_mip->w2);

        while (!xmen_mip->ERNUM && EasyXMen_sm2_size(xmen_mip->w2) != 0)
        { /* use binary method */
            if (XMen_Compare(xmen_mip->w1, xmen_mip->w2) >= 0)
            { /* vp=v*vp-p, v=v*v-2 */
                nres_modmult(xmen_mip->w8, xmen_mip->w9, xmen_mip->w8);
                nres_modsub(xmen_mip->w8, xmen_mip->w3, xmen_mip->w8);
                nres_modmult(xmen_mip->w9, xmen_mip->w9, xmen_mip->w9);
                nres_modsub(xmen_mip->w9, xmen_mip->w4, xmen_mip->w9);
                XMen_Subtract(xmen_mip->w1, xmen_mip->w2, xmen_mip->w1);
            }
            else
            { /* v=v*vp-p, vp=vp*vp-2 */
                nres_modmult(xmen_mip->w9, xmen_mip->w8, xmen_mip->w9);
                nres_modsub(xmen_mip->w9, xmen_mip->w3, xmen_mip->w9);
                nres_modmult(xmen_mip->w8, xmen_mip->w8, xmen_mip->w8);
                nres_modsub(xmen_mip->w8, xmen_mip->w4, xmen_mip->w8);
            }
            subdiv(xmen_mip->w2, 2, xmen_mip->w2);
        }
    }
#endif

    copy(xmen_mip->w9, v);
    if (v != vp)
        copy(xmen_mip->w8, vp);
    XMEN_OUT
}

/**
 * @description Performs modular division using n-residues.
 *        This function calculates $ w = x / y \mod n $ and stores the result in `w`.
 *        It handles the case where the division is not possible by checking the greatest common divisor (GCD).
 * @param[in] x   : Pointer to the multi-precision integer representing the dividend.
 * @param[in] y   : Pointer to the multi-precision integer representing the divisor.
 * @param[out] w  : Pointer to the multi-precision integer where the result of $ x / y \mod n $ will be stored.
 * @return int    : Returns the GCD of `y` and `n`. If the GCD is not 1, it indicates that the division failed.
 * @retval 0      : Division failed or an error occurred.
 */
static int nres_moddiv(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr w)
{ /* Modular division using n-residues w=x/y mod n */
    int gcd;
    if (xmen_mip->ERNUM)
        return 0;

    XMEN_IN(85)

    if (x == y)
    { /* Illegal parameter usage */
        EasyXMen_berror(XMEN_ERR_BAD_PARAMETERS);
        XMEN_OUT

        return 0;
    }
    redc(y, xmen_mip->w6);
    gcd = invmodp(xmen_mip->w6, xmen_mip->modulus, xmen_mip->w6);

    if (gcd != 1)
        zero(w); /* fails silently and returns 0 */
    else
    {
        nres(xmen_mip->w6, xmen_mip->w6);
        nres_modmult(x, xmen_mip->w6, w);
    }
    XMEN_OUT
    return gcd;
}

/**
 * @description Computes the modular square root of a multi-precision integer.
 *        This function calculates $ w = \sqrt{x} \mod p $, where `p` is assumed to be prime.
 *        It handles special cases and uses different algorithms based on the properties of the modulus.
 * @param[in] x   : Pointer to the multi-precision integer whose square root is to be calculated.
 * @param[out] w  : Pointer to the multi-precision integer where the result (square root) will be stored.
 * @return BOOL   : Returns TRUE if the square root was successfully computed, otherwise returns FALSE.
 * @retval TRUE   : The square root was successfully computed.
 * @retval FALSE  : The square root could not be computed or an error occurred.
 */
static BOOL nres_sqroot(XMen_BigTPtr x, XMen_BigTPtr w)
{ /* w=sqrt(x) mod p. This depends on p being prime! */
    int t, js;
    if (xmen_mip->ERNUM)
        return FALSE;

    copy(x, w);
    if (EasyXMen_sm2_size(w) == 0)
        return TRUE;

    XMEN_IN(100)

    redc(w, w); /* get it back into normal form */

    if (EasyXMen_sm2_size(w) == 1) /* square root of 1 is 1 */
    {
        nres(w, w);
        XMEN_OUT
        return TRUE;
    }

    if (EasyXMen_sm2_size(w) == 4) /* square root of 4 is 2 */
    {
        convert(2, w);
        nres(w, w);
        XMEN_OUT
        return TRUE;
    }

    if (jack(w, xmen_mip->modulus) != 1)
    { /* Jacobi test */
        zero(w);
        XMEN_OUT
        return FALSE;
    }

    js = xmen_mip->pmod8 % 4 - 2; /* 1 mod 4 or 3 mod 4 prime? */

    XMen_Incr(xmen_mip->modulus, js, xmen_mip->w10);
    subdiv(xmen_mip->w10, 4, xmen_mip->w10); /* (p+/-1)/4 */

    if (js == 1)
    { /* 3 mod 4 primes - do a quick and dirty sqrt(x)=x^(p+1)/4 mod p */
        nres(w, xmen_mip->w2);
        copy(xmen_mip->one, w);
        for (;;)
        { /* Simple Right-to-Left exponentiation */

            if (xmen_mip->user != NULL)
                (*xmen_mip->user)();
            if (subdiv(xmen_mip->w10, 2, xmen_mip->w10) != 0)
                nres_modmult(w, xmen_mip->w2, w);
            if (xmen_mip->ERNUM || EasyXMen_sm2_size(xmen_mip->w10) == 0)
                break;
            nres_modmult(xmen_mip->w2, xmen_mip->w2, xmen_mip->w2);
        }

        /*     nres_moddiv(xmen_mip->one,w,xmen_mip->w11);
               nres_modadd(xmen_mip->w11,w,xmen_mip->w3);
               nres_lucas(xmen_mip->w3,xmen_mip->w10,w,w);
               nres_modadd(xmen_mip->w11,xmen_mip->one,xmen_mip->w11);
               nres_moddiv(w,xmen_mip->w11,w); */
    }
    else
    { /* 1 mod 4 primes */
        for (t = 1;; t++)
        { /* t=1.5 on average */
            if (t == 1)
                copy(w, xmen_mip->w4);
            else
            {
                premult(w, t, xmen_mip->w4);
                XMen_Divide(xmen_mip->w4, xmen_mip->modulus, xmen_mip->modulus);
                premult(xmen_mip->w4, t, xmen_mip->w4);
                XMen_Divide(xmen_mip->w4, xmen_mip->modulus, xmen_mip->modulus);
            }

            decr(xmen_mip->w4, 4, xmen_mip->w1);
            if (jack(xmen_mip->w1, xmen_mip->modulus) == js)
                break;
            if (xmen_mip->ERNUM)
                break;
        }

        decr(xmen_mip->w4, 2, xmen_mip->w3);
        nres(xmen_mip->w3, xmen_mip->w3);
        nres_lucas(xmen_mip->w3, xmen_mip->w10, w, w); /* heavy lifting done here */
        if (t != 1)
        {
            convert(t, xmen_mip->w11);
            nres(xmen_mip->w11, xmen_mip->w11);
            nres_moddiv(w, xmen_mip->w11, w);
        }
    }

    XMEN_OUT
    return TRUE;
}

/**
 * @description Generates a random number using the Marsaglia & Zaman algorithm.
 *        This function produces a pseudo-random number by manipulating an internal state array.
 *        The size of the underlying type determines how the state is accessed and combined to form the output.
 * @return XMen_SmallType : The generated pseudo-random number.
 */
static XMen_SmallType brand(void)
{ /* Marsaglia & Zaman random number generator */
    int i, k;
    unsigned int pdiff, t;
    XMen_SmallType r;
    if (xmen_mip->lg2b > 32)
    { /* underlying type is > 32 bits. Assume <= 64 bits */
        xmen_mip->rndptr += 2;
        if (xmen_mip->rndptr < NK - 1)
        {
            r = (XMen_SmallType)xmen_mip->ira[xmen_mip->rndptr];
            r = xmen_shiftbits(r, xmen_mip->lg2b - 32);
            r += (XMen_SmallType)xmen_mip->ira[xmen_mip->rndptr + 1];
            return r;
        }
    }
    else
    {
        xmen_mip->rndptr++;
        if (xmen_mip->rndptr < NK)
            return (XMen_SmallType)xmen_mip->ira[xmen_mip->rndptr];
    }
    xmen_mip->rndptr = 0;
    for (i = 0, k = NK - NJ; i < NK; i++, k++)
    { /* calculate next NK values */
        if (k == NK)
            k = 0;
        t = xmen_mip->ira[k];
        pdiff = t - xmen_mip->ira[i] - xmen_mip->borrow;
        if (pdiff < t)
            xmen_mip->borrow = 0;
        if (pdiff > t)
            xmen_mip->borrow = 1;
        xmen_mip->ira[i] = pdiff;
    }
    if (xmen_mip->lg2b > 32)
    { /* double up */
        r = (XMen_SmallType)xmen_mip->ira[0];
        r = xmen_shiftbits(r, xmen_mip->lg2b - 32);
        r += (XMen_SmallType)xmen_mip->ira[1];
        return r;
    }
    else
        return (XMen_SmallType)(xmen_mip->ira[0]);
}

/**
 * @description Copies the contents of one elliptic curve point to another.
 *        This function copies the X, Y coordinates and optionally the Z coordinate if the point is in general form.
 *        It also copies the marker indicating the type or state of the point.
 * @param[in] a   : Pointer to the source elliptic curve point.
 * @param[out] b  : Pointer to the destination elliptic curve point.
 * @return void
 */
static void epoint_copy(XMen_ECurvePointT *a, XMen_ECurvePointT *b)
{
    if (a == b || b == NULL)
        return;

    copy(a->X, b->X);
    copy(a->Y, b->Y);
    if (a->marker == XMEN_EPOINT_GENERAL)
        copy(a->Z, b->Z);
    b->marker = a->marker;
    return;
}

/**
 * @description Initializes an elliptic curve point using pre-allocated memory.
 *        This function sets up an elliptic curve point structure within a specified memory block,
 *        initializing its components based on the coordinate system used.
 * @param[in] mem : Pointer to the pre-allocated memory block.
 * @param[in] index : Index used to calculate the offset within the memory block.
 * @param[in] sz : Size parameter for memory allocation.
 * @return XMen_ECurvePointT* : Pointer to the initialized elliptic curve point.
 */
static XMen_ECurvePointT *epoint_init_mem_variable(char *mem, int index, int sz)
{
    XMen_ECurvePointT *p;
    char *ptr;
    int offset, r;

    offset = 0;
    r = (unsigned long)mem % XMEN_SL;
    if (r > 0)
        offset = XMEN_SL - r;

    if (xmen_mip->coord == XMEN_AFFINE)
        p = (XMen_ECurvePointT *)&mem[offset + index * xmen_esize_a(sz)];
    else
        p = (XMen_ECurvePointT *)&mem[offset + index * xmen_esize(sz)];

    ptr = (char *)p + sizeof(XMen_ECurvePointT);
    p->X = mirvar_mem_variable(ptr, 0, sz);
    p->Y = mirvar_mem_variable(ptr, 1, sz);
    if (xmen_mip->coord != XMEN_AFFINE)
        p->Z = mirvar_mem_variable(ptr, 2, sz);
    p->marker = XMEN_EPOINT_INFINITY;
    return p;
}

/**
 * @description Negates an elliptic curve point.
 *        This function negates the Y-coordinate of the point, effectively reflecting it over the x-axis.
 *        Points at infinity remain unchanged.
 * @param[in,out] p : Pointer to the elliptic curve point to negate.
 * @return void
 */
static void epoint_negate(XMen_ECurvePointT *p)
{ /* negate a point */
    if (xmen_mip->ERNUM)
        return;
    if (p->marker == XMEN_EPOINT_INFINITY)
        return;

    XMEN_IN(121)
    if (EasyXMen_sm2_size(p->Y) != 0)
        xmen_psub(xmen_mip->modulus, p->Y, p->Y);
    XMEN_OUT
}

/**
 * @description Initializes an elliptic curve point using pre-allocated memory.
 *        This function initializes an elliptic curve point within a specified memory block,
 *        utilizing the `epoint_init_mem_variable` helper function with a size parameter derived from the context.
 * @param[in] mem : Pointer to the pre-allocated memory block.
 * @param[in] index : Index used to calculate the offset within the memory block.
 * @return XMen_ECurvePointT* : Pointer to the initialized elliptic curve point, or NULL if an error occurred.
 * @retval NULL : Initialization failed due to an error condition.
 */
static XMen_ECurvePointT *epoint_init_mem(char *mem, int index)
{
    if (xmen_mip->ERNUM)
        return NULL;

    return epoint_init_mem_variable(mem, index, xmen_mip->nib - 1);
}
#ifndef XMEN_STATIC

/* allocate space for a number of epoints from the heap */
/**
 * @description Allocates memory for elliptic curve point operations.
 *        This function reserves memory for a specified number of elliptic curve points,
 *        taking into account the coordinate system (affine or projective) used by the context.
 * @param[in] num : The number of elliptic curve points to allocate memory for.
 * @return void* : Pointer to the allocated memory block, or NULL if allocation fails.
 */
static void *ecp_memalloc(int num)
{

    if (xmen_mip->coord == XMEN_AFFINE)
        return EasyXMen_alloc(xmen_ecp_reserve_a(num, xmen_mip->nib - 1), 1);
    else
        return EasyXMen_alloc(xmen_ecp_reserve(num, xmen_mip->nib - 1), 1);
}

#endif
/**
 * @description Computes the modular negation of a multi-precision integer.
 *        This function calculates $ w = -x \mod n $ and stores the result in `w`.
 *        It handles special cases, such as when `x` is zero, and uses optimized routines if available.
 * @param[in] x   : Pointer to the multi-precision integer to negate.
 * @param[out] w  : Pointer to the multi-precision integer where the negated value will be stored.
 * @return void
 */
static void nres_negate(XMen_BigTPtr x, XMen_BigTPtr w)
{
    if (EasyXMen_sm2_size(x) == 0)
    {
        zero(w);
        return;
    }
#ifdef XMEN_COMBA
    if (xmen_mip->ACTIVE)
    {
        comba_negate(x, w);
        return;
    }
    else
    {
#endif
        if (xmen_mip->ERNUM)
            return;

        XMEN_IN(92)
        xmen_psub(xmen_mip->modulus, x, w);
        XMEN_OUT

#ifdef XMEN_COMBA
    }
#endif
}

/**
 * @description Performs a modular division by 2 on a multi-precision integer.
 *        This function calculates $ w = x / 2 \mod n $ and handles cases where `x` is odd by adjusting it using the modulus.
 * @param[in] x   : Pointer to the multi-precision integer to divide by 2.
 * @param[out] w  : Pointer to the multi-precision integer where the result will be stored.
 * @return void
 */
static void nres_div2(XMen_BigTPtr x, XMen_BigTPtr w)
{

    XMEN_IN(198)
    copy(x, xmen_mip->w1);
    if (remain(xmen_mip->w1, 2) != 0)
        XMen_Add(xmen_mip->w1, xmen_mip->modulus, xmen_mip->w1);
    subdiv(xmen_mip->w1, 2, xmen_mip->w1);
    copy(xmen_mip->w1, w);

    XMEN_OUT
}

/**
 * @description Adds two elliptic curve points on the active elliptic curve.
 *        This function performs point addition $ pa += p $ on an elliptic curve,
 *        handling both affine and projective coordinate systems. It ensures correct operation even when one of the points is not normalized.
 * @param[in] p   : Pointer to the first elliptic curve point to add.
 * @param[in,out] pa : Pointer to the second elliptic curve point, which will be updated with the sum.
 * @return BOOL   : Returns TRUE if the addition was successful, otherwise returns FALSE.
 * @retval TRUE   : The addition was successful.
 * @retval FALSE  : An error occurred during addition or the operation should have been a doubling.
 */
static BOOL ecurve_padd(XMen_ECurvePointT *p, XMen_ECurvePointT *pa)
{ /* primitive add two epoints on the active ecurve - pa+=p;   *
   * note that if p is normalized, its Z coordinate isn't used */

    if (xmen_mip->coord == XMEN_AFFINE)
    { /* 1 sqr, 1 mul, 1 div */
        nres_modsub(p->Y, pa->Y, xmen_mip->w8);
        nres_modsub(p->X, pa->X, xmen_mip->w6);
        if (EasyXMen_sm2_size(xmen_mip->w6) == 0)
        { /* divide by 0 */
            if (EasyXMen_sm2_size(xmen_mip->w8) == 0)
            { /* should have doubled ! */
                return FALSE;
            }
            else
            { /* point at infinity */
                XMen_EpointSet(NULL, NULL, 0, pa);
                return TRUE;
            }
        }
        if (nres_moddiv(xmen_mip->w8, xmen_mip->w6, xmen_mip->w8) > 1)
        {
            XMen_EpointSet(NULL, NULL, 0, pa);
            EasyXMen_berror(XMEN_ERR_COMPOSITE_MODULUS);
            return TRUE;
        }

        nres_modmult(xmen_mip->w8, xmen_mip->w8, xmen_mip->w2); /* w2=m^2 */
        nres_modsub(xmen_mip->w2, p->X, xmen_mip->w1);          /* w1=m^2-x1-x2 */
        nres_modsub(xmen_mip->w1, pa->X, xmen_mip->w1);

        nres_modsub(pa->X, xmen_mip->w1, xmen_mip->w2);
        nres_modmult(xmen_mip->w2, xmen_mip->w8, xmen_mip->w2);
        nres_modsub(xmen_mip->w2, pa->Y, pa->Y);
        copy(xmen_mip->w1, pa->X);

        pa->marker = XMEN_EPOINT_NORMALIZED;
        return TRUE;
    }

    if (p->marker != XMEN_EPOINT_NORMALIZED)
    {
        nres_modmult(p->Z, p->Z, xmen_mip->w6);
        nres_modmult(pa->X, xmen_mip->w6, xmen_mip->w1);
        nres_modmult(xmen_mip->w6, p->Z, xmen_mip->w6);
        nres_modmult(pa->Y, xmen_mip->w6, xmen_mip->w8);
    }
    else
    {
        copy(pa->X, xmen_mip->w1);
        copy(pa->Y, xmen_mip->w8);
    }
    if (pa->marker == XMEN_EPOINT_NORMALIZED)
        copy(xmen_mip->one, xmen_mip->w6);

    else
        nres_modmult(pa->Z, pa->Z, xmen_mip->w6);
    nres_modmult(p->X, xmen_mip->w6, xmen_mip->w4);
    if (pa->marker != XMEN_EPOINT_NORMALIZED)
        nres_modmult(xmen_mip->w6, pa->Z, xmen_mip->w6);
    nres_modmult(p->Y, xmen_mip->w6, xmen_mip->w5);
    nres_modsub(xmen_mip->w1, xmen_mip->w4, xmen_mip->w1);
    nres_modsub(xmen_mip->w8, xmen_mip->w5, xmen_mip->w8);

    /* w8 contains the numerator of the slope */

    if (EasyXMen_sm2_size(xmen_mip->w1) == 0)
    {
        if (EasyXMen_sm2_size(xmen_mip->w8) == 0)
        { /* should have doubled ! */
            return FALSE;
        }
        else
        { /* point at infinity */
            XMen_EpointSet(NULL, NULL, 0, pa);
            return TRUE;
        }
    }
    nres_modadd(xmen_mip->w4, xmen_mip->w4, xmen_mip->w6);
    nres_modadd(xmen_mip->w1, xmen_mip->w6, xmen_mip->w4);
    nres_modadd(xmen_mip->w5, xmen_mip->w5, xmen_mip->w6);
    nres_modadd(xmen_mip->w8, xmen_mip->w6, xmen_mip->w5);

    if (p->marker != XMEN_EPOINT_NORMALIZED)
    {
        if (pa->marker != XMEN_EPOINT_NORMALIZED)
            nres_modmult(pa->Z, p->Z, xmen_mip->w3);
        else
            copy(p->Z, xmen_mip->w3);
        nres_modmult(xmen_mip->w3, xmen_mip->w1, pa->Z);
    }
    else
    {
        if (pa->marker != XMEN_EPOINT_NORMALIZED)
            nres_modmult(pa->Z, xmen_mip->w1, pa->Z);
        else
            copy(xmen_mip->w1, pa->Z);
    }
    nres_modmult(xmen_mip->w1, xmen_mip->w1, xmen_mip->w6);
    nres_modmult(xmen_mip->w1, xmen_mip->w6, xmen_mip->w1);
    nres_modmult(xmen_mip->w6, xmen_mip->w4, xmen_mip->w6);
    nres_modmult(xmen_mip->w8, xmen_mip->w8, xmen_mip->w4);

    nres_modsub(xmen_mip->w4, xmen_mip->w6, pa->X);
    nres_modsub(xmen_mip->w6, pa->X, xmen_mip->w6);
    nres_modsub(xmen_mip->w6, pa->X, xmen_mip->w6);
    nres_modmult(xmen_mip->w8, xmen_mip->w6, xmen_mip->w2);
    nres_modmult(xmen_mip->w1, xmen_mip->w5, xmen_mip->w1);
    nres_modsub(xmen_mip->w2, xmen_mip->w1, xmen_mip->w5);

    /* divide by 2 */

    nres_div2(xmen_mip->w5, pa->Y);

    pa->marker = XMEN_EPOINT_GENERAL;
    return TRUE;
}

/**
 * @description Adds two elliptic curve points on the active elliptic curve.
 *        This function performs point addition $ pa = pa + p $ on an elliptic curve,
 *        handling special cases such as point doubling and adding the point at infinity.
 *        It returns a status code indicating the type of operation performed or if an error occurred.
 * @param[in] p   : Pointer to the first elliptic curve point to add.
 * @param[in,out] pa : Pointer to the second elliptic curve point, which will be updated with the sum.
 * @return int    : Status code indicating the result of the operation.
 * @retval XMEN_ADD : The addition was successful.
 * @retval XMEN_DOUBLE : A point doubling operation was performed instead of addition.
 * @retval XMEN_OVER : An error occurred, such as encountering a point at infinity or another exceptional condition.
 */
static int ecurve_add(XMen_ECurvePointT *p, XMen_ECurvePointT *pa)
{ /* pa=pa+p; */
    /* An ephemeral pointer to the line slope is returned */

    if (xmen_mip->ERNUM)
        return XMEN_OVER;

    XMEN_IN(94)

    if (p == pa)
    {
        ecurve_double(pa);
        XMEN_OUT
        if (pa->marker == XMEN_EPOINT_INFINITY)
            return XMEN_OVER;
        return XMEN_DOUBLE;
    }
    if (pa->marker == XMEN_EPOINT_INFINITY)
    {
        epoint_copy(p, pa);
        XMEN_OUT
        return XMEN_ADD;
    }
    if (p->marker == XMEN_EPOINT_INFINITY)
    {
        XMEN_OUT
        return XMEN_ADD;
    }

    if (!ecurve_padd(p, pa))
    {
        ecurve_double(pa);
        XMEN_OUT
        return XMEN_DOUBLE;
    }
    XMEN_OUT
    if (pa->marker == XMEN_EPOINT_INFINITY)
        return XMEN_OVER;
    return XMEN_ADD;
}

/**
 * @description Computes the modular inverses of multiple multi-precision integers.
 *        This function calculates $ w[i] = 1 / x[i] \mod n $ for a set of multi-precision integers,
 *        ensuring that the input and output arrays are distinct.
 * @param[in] m   : The number of multi-precision integers to compute inverses for.
 * @param[in] x   : Pointer to an array of multi-precision integers for which to find inverses.
 * @param[out] w  : Pointer to an array where the computed inverses will be stored.
 * @return BOOL   : Returns TRUE if the operation was successful, otherwise returns FALSE.
 * @retval TRUE   : All inverses were successfully computed.
 * @retval FALSE  : An error occurred, such as invalid parameters or division by zero.
 */
static BOOL nres_multi_inverse(int m, XMen_BigTPtr *x, XMen_BigTPtr *w)
{ /* find w[i]=1/x[i] mod n, for i=0 to m-1 *
   * x and w MUST be distinct               */
    int i;
    if (m == 0)
        return TRUE;
    if (m < 0)
        return FALSE;
    XMEN_IN(118)

    if (x == w)
    {
        EasyXMen_berror(XMEN_ERR_BAD_PARAMETERS);
        XMEN_OUT
        return FALSE;
    }

    if (m == 1)
    {
        copy(xmen_mip->one, w[0]);
        nres_moddiv(w[0], x[0], w[0]);
        XMEN_OUT
        return TRUE;
    }

    convert(1, w[0]);
    copy(x[0], w[1]);
    for (i = 2; i < m; i++)
        nres_modmult(w[i - 1], x[i - 1], w[i]);

    nres_modmult(w[m - 1], x[m - 1], xmen_mip->w6); /* y=x[0]*x[1]*x[2]....x[m-1] */
    if (EasyXMen_sm2_size(xmen_mip->w6) == 0)
    {
        EasyXMen_berror(XMEN_ERR_DIV_BY_ZERO);
        XMEN_OUT
        return FALSE;
    }

    redc(xmen_mip->w6, xmen_mip->w6);
    redc(xmen_mip->w6, xmen_mip->w6);

    invmodp(xmen_mip->w6, xmen_mip->modulus, xmen_mip->w6);

    /* Now y=1/y */

    copy(x[m - 1], xmen_mip->w5);
    nres_modmult(w[m - 1], xmen_mip->w6, w[m - 1]);

    for (i = m - 2;; i--)
    {
        if (i == 0)
        {
            nres_modmult(xmen_mip->w5, xmen_mip->w6, w[0]);
            break;
        }
        nres_modmult(w[i], xmen_mip->w5, w[i]);
        nres_modmult(w[i], xmen_mip->w6, w[i]);
        nres_modmult(xmen_mip->w5, x[i], xmen_mip->w5);
    }

    XMEN_OUT
    return TRUE;
}

/**
 * @description Normalizes an array of elliptic curve points.
 *        This function normalizes a set of elliptic curve points, converting them from projective to affine coordinates,
 *        using a workspace array for intermediate calculations. It ensures that the normalization process is only applied
 *        to points in projective coordinates and handles up to a maximum number of points defined by `XMEN_MAX_M_T_S`.
 * @param[in] m   : The number of elliptic curve points to normalize.
 * @param[in,out] work : Pointer to an array used as workspace for intermediate calculations.
 * @param[in,out] p    : Pointer to an array of pointers to elliptic curve points to be normalized.
 * @return BOOL   : Returns TRUE if the operation was successful, otherwise returns FALSE.
 * @retval TRUE   : All points were successfully normalized.
 * @retval FALSE  : An error occurred, such as invalid parameters or exceeding the maximum number of points.
 */
static BOOL epoint_multi_norm(int m, XMen_BigTPtr *work, XMen_ECurvePointT **p)
{ /* Normalise an array of points of length m<XMEN_MAX_M_T_S - requires a workspace array of length m */

    int i;
    XMen_BigTPtr w[XMEN_MAX_M_T_S];
    if (xmen_mip->coord == XMEN_AFFINE)
        return TRUE;
    if (xmen_mip->ERNUM)
        return FALSE;
    if (m > XMEN_MAX_M_T_S)
        return FALSE;

    XMEN_IN(190)

    for (i = 0; i < m; i++)
    {
        if (p[i]->marker == XMEN_EPOINT_NORMALIZED)
            w[i] = xmen_mip->one;
        else
            w[i] = p[i]->Z;
    }

    if (!nres_multi_inverse(m, w, work))
    {
        XMEN_OUT
        return FALSE;
    }

    for (i = 0; i < m; i++)
    {
        copy(xmen_mip->one, p[i]->Z);
        p[i]->marker = XMEN_EPOINT_NORMALIZED;
        nres_modmult(work[i], work[i], xmen_mip->w1);
        nres_modmult(p[i]->X, xmen_mip->w1, p[i]->X); /* X/ZZ */
        nres_modmult(xmen_mip->w1, work[i], xmen_mip->w1);
        nres_modmult(p[i]->Y, xmen_mip->w1, p[i]->Y); /* Y/ZZZ */
    }
    XMEN_OUT
    return TRUE;
}

/**
 * @description Computes a sliding window Non-Adjacent Form (NAF) value for multi-precision integers.
 *        This function calculates the NAF value using fractional windows, where "store" precomputed values are used.
 *        It starts scanning from the i-th bit of `x` and combines `x` and `x3` (which is 3*x) to produce the NAF.
 *        The function returns the sliding window value along with the number of bits processed (`nbs`) and the number
 *        of trailing zeros detected (`nzs`). The return value is zero or odd, within the range -21 <= n <= +21.
 * @param[in] x   : Pointer to the multi-precision integer to process.
 * @param[in] x3  : Pointer to the multi-precision integer representing 3*x.
 * @param[in] i   : The starting bit position for scanning.
 * @param[out] nbs: Pointer to an integer where the number of bits processed will be stored.
 * @param[out] nzs: Pointer to an integer where the number of trailing zeros detected will be stored.
 * @param[in] store : The number of precomputed values available for use in the calculation.
 * @return int    : Returns the sliding window NAF value.
 * @retval n      : The computed NAF value, which is zero or odd and within the range -21 <= n <= +21.
 */
static int xmen_naf_window(XMen_BigTPtr x, XMen_BigTPtr x3, int i, int *nbs, int *nzs, int store)
{ /* returns sliding window value, using fractional windows   *
   * where "store" precomputed values are precalulated and    *
   * stored. Scanning starts at the i-th bit of  x. nbs is    *
   * the number of bits processed. nzs is number of           *
   * additional trailing zeros detected. x and x3 (which is   *
   * 3*x) are combined to produce the NAF (non-adjacent       *
   * form). So if x=11011(27) and x3 is 1010001, the LSB is   *
   * ignored and the value 100T0T (32-4-1=27) processed,      *
   * where T is -1. Note x.P = (3x-x)/2.P. This value will    *
   * return +7, with nbs=4 and nzs=1, having stopped after    *
   * the first 4 bits. If it goes too far, it must backtrack  *
   * Note in an NAF non-zero elements are never side by side, *
   * so 10T10T won't happen. NOTE: return value n zero or     *
   * odd, -21 <= n <= +21     */

    int nb, j, r, biggest;

    /* get first bit */
    nb = xmen_testbit(x3, i) - xmen_testbit(x, i);

    *nbs = 1;
    *nzs = 0;
    if (nb == 0)
        return 0;
    if (i == 0)
        return nb;

    biggest = 2 * store - 1;

    if (nb > 0)
        r = 1;
    else
        r = (-1);

    for (j = i - 1; j > 0; j--)
    {
        (*nbs)++;
        r *= 2;
        nb = xmen_testbit(x3, j) - xmen_testbit(x, j);
        if (nb > 0)
            r += 1;
        if (nb < 0)
            r -= 1;
        if (abs(r) > biggest)
            break;
    }

    if (r % 2 != 0 && j != 0)
    { /* backtrack */
        if (nb > 0)
            r = (r - 1) / 2;
        if (nb < 0)
            r = (r + 1) / 2;
        (*nbs)--;
    }

    while (r % 2 == 0)
    { /* remove trailing zeros */
        r /= 2;
        (*nzs)++;
        (*nbs)--;
    }
    return r;
}

/**
 * @description Subtracts one elliptic curve point from another on the active elliptic curve.
 *        This function performs point subtraction $ pa = pa - p $ on an elliptic curve,
 *        handling special cases such as subtracting the same point or a point at infinity.
 *        It leverages point negation and addition to achieve subtraction.
 * @param[in] p   : Pointer to the elliptic curve point to subtract.
 * @param[in,out] pa : Pointer to the elliptic curve point from which `p` will be subtracted.
 * @return int    : Status code indicating the result of the operation.
 * @retval XMEN_ADD : The subtraction was successful, equivalent to performing an addition.
 * @retval XMEN_OVER : An error occurred, such as encountering a point at infinity or another exceptional condition.
 */
static int ecurve_sub(XMen_ECurvePointT *p, XMen_ECurvePointT *pa)
{
    int r;
    if (xmen_mip->ERNUM)
        return XMEN_OVER;

    XMEN_IN(104)

    if (p == pa)
    {
        XMen_EpointSet(NULL, NULL, 0, pa);
        XMEN_OUT
        return XMEN_OVER;
    }
    if (p->marker == XMEN_EPOINT_INFINITY)
    {
        XMEN_OUT
        return XMEN_ADD;
    }

    epoint_negate(p);
    r = ecurve_add(p, pa);
    epoint_negate(p);

    XMEN_OUT
    return r;
}

/**
 * @description Securely clears and optionally frees a block of memory used for elliptic curve points.
 *        This function securely wipes the memory area containing elliptic curve point data,
 *        ensuring that sensitive information is not left in memory. It also frees the memory if
 *        dynamic allocation is enabled.
 * @param[in] mem : Pointer to the memory block to clear and free.
 * @param[in] num : The number of elements (elliptic curve points) in the memory block.
 * @return void
 */
static void ecp_memkill(char *mem, int num)
{
    if (mem == NULL)
        return;

    if (xmen_mip->coord == XMEN_AFFINE)
        memset(mem, 0, xmen_ecp_reserve_a(num, xmen_mip->nib - 1));
    else
        memset(mem, 0, xmen_ecp_reserve(num, xmen_mip->nib - 1));

#ifndef XMEN_STATIC
    EasyXMen_free(mem);
#endif
}

/**
 * @description Securely clears and optionally frees a block of memory.
 *        This function securely wipes the memory area, ensuring that sensitive information is not left in memory.
 *        It also frees the memory if dynamic allocation is enabled.
 * @param[in] mem : Pointer to the memory block to clear and free.
 * @param[in] len : The length of the memory block in bytes.
 * @return void
 */
static void memkill(char *mem, int len)
{
    if (mem == NULL)
        return;
    memset(mem, 0, xmen_big_reserve(len, xmen_mip->nib - 1));
#ifndef XMEN_STATIC
    EasyXMen_free(mem);
#endif
}

/**
 * @description Normalizes an elliptic curve point from projective to affine coordinates.
 *        This function converts a point from projective to affine coordinates by dividing the X and Y components
 *        by the Z component's appropriate power. If normalization fails, it sets the point to the point at infinity.
 * @param[in,out] p : Pointer to the elliptic curve point to normalize.
 * @return BOOL     : Returns TRUE if the normalization was successful, otherwise returns FALSE.
 * @retval TRUE     : The point was successfully normalized.
 * @retval FALSE    : An error occurred during normalization.
 */
static BOOL epoint_norm(XMen_ECurvePointT *p)
{ /* normalise a point */

    if (xmen_mip->coord == XMEN_AFFINE)
        return TRUE;
    if (p->marker != XMEN_EPOINT_GENERAL)
        return TRUE;

    if (xmen_mip->ERNUM)
        return FALSE;

    XMEN_IN(117)

    copy(xmen_mip->one, xmen_mip->w8);

    if (nres_moddiv(xmen_mip->w8, p->Z, xmen_mip->w8) > 1) /* 1/Z  */
    {
        XMen_EpointSet(NULL, NULL, 0, p);
        EasyXMen_berror(XMEN_ERR_COMPOSITE_MODULUS);
        XMEN_OUT
        return FALSE;
    }

    nres_modmult(xmen_mip->w8, xmen_mip->w8, xmen_mip->w1); /* 1/ZZ */
    nres_modmult(p->X, xmen_mip->w1, p->X);                 /* X/ZZ */
    nres_modmult(xmen_mip->w1, xmen_mip->w8, xmen_mip->w1); /* 1/ZZZ */
    nres_modmult(p->Y, xmen_mip->w1, p->Y);                 /* Y/ZZZ */

    copy(xmen_mip->one, p->Z);

    p->marker = XMEN_EPOINT_NORMALIZED;
    XMEN_OUT

    return TRUE;
}

/**
 * @description Sets the size of the I/O buffer and reallocates it accordingly.
 *        This function resizes the I/O buffer to the specified length, initializing it with zeros,
 *        and ensures that any existing buffer is properly freed before allocating a new one.
 *        If the specified length is zero, the current buffer is freed and no new buffer is allocated.
 * @param[in] len : The desired size of the I/O buffer in bytes. Must be non-negative.
 * @return void
 */
static void set_io_buffer_size(int len)
{
    int i;
    if (len < 0)
        return;
    XMEN_IN(142)
    for (i = 0; i < xmen_mip->IOBSIZ; i++)
        xmen_mip->IOBUFF[i] = 0;
    EasyXMen_free(xmen_mip->IOBUFF);
    if (len == 0)
    {
        XMEN_OUT
        return;
    }
    xmen_mip->IOBSIZ = len;
    xmen_mip->IOBUFF = (char *)EasyXMen_alloc(len + 1, 1);
    xmen_mip->IOBUFF[0] = '\0';
    XMEN_OUT
}

//*Solinas's Joint Sparse Form */
/**
 * @description Adds $2^n$ to a multi-precision integer assuming the nth bit is zero.
 *        This function sets the nth bit of the multi-precision integer `x`, effectively adding $2^n$ to `x`.
 *        It should be used with caution as it assumes the nth bit is initially zero.
 * @param[in,out] x : Pointer to the multi-precision integer to modify.
 * @param[in] n     : The position of the bit to set (the exponent for $2^n$).
 * @return void
 * @note Use with care! It is assumed that the nth bit of `x` is initially zero.
 */
static void xmen_addbit(XMen_BigTPtr x, int n)
{ /* add 2^n to positive x - where you know that bit is zero. Use with care! */
    unsigned int m = n / xmen_mip->lg2b;
    x->w[m] += xmen_shiftbits((XMen_SmallType)1, n % xmen_mip->lg2b);
    if (x->len < m + 1)
        x->len = m + 1;
}

/**
 * @description Computes the Joint Sparse Form (JSF) for two multi-precision integers.
 *        This function calculates the JSF representation of two multi-precision integers `k0` and `k1`,
 *        storing the results in `u0p`, `u0m`, `u1p`, and `u1m`. The JSF is a method to represent two integers
 *        in a way that minimizes the number of non-zero digits, which can be useful for optimizing certain
 *        cryptographic operations. The function handles both binary and non-binary bases.
 * @param[in] k0    : Pointer to the first multi-precision integer.
 * @param[in] k1    : Pointer to the second multi-precision integer.
 * @param[out] u0p  : Pointer to the multi-precision integer where positive parts of the JSF for `k0` will be stored.
 * @param[out] u0m  : Pointer to the multi-precision integer where negative parts of the JSF for `k0` will be stored.
 * @param[out] u1p  : Pointer to the multi-precision integer where positive parts of the JSF for `k1` will be stored.
 * @param[out] u1m  : Pointer to the multi-precision integer where negative parts of the JSF for `k1` will be stored.
 * @return void
 */
static void xmen_jsf(XMen_BigTPtr k0, XMen_BigTPtr k1, XMen_BigTPtr u0p, XMen_BigTPtr u0m, XMen_BigTPtr u1p, XMen_BigTPtr u1m)
{
    int j, u0, u1, d0, d1, l0, l1;
    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(191)

    d0 = d1 = 0;

    convert(1, xmen_mip->w1);
    copy(k0, xmen_mip->w2);
    copy(k1, xmen_mip->w3);
    zero(u0p);
    zero(u0m);
    zero(u1p);
    zero(u1m);

    j = 0;
    while (!xmen_mip->ERNUM)
    {
        if (EasyXMen_sm2_size(xmen_mip->w2) == 0 && d0 == 0 && EasyXMen_sm2_size(xmen_mip->w3) == 0 && d1 == 0)
            break;
        l0 = remain(xmen_mip->w2, 8);
        l0 = (l0 + d0) & 0x7;
        l1 = remain(xmen_mip->w3, 8);
        l1 = (l1 + d1) & 0x7;

        if (l0 % 2 == 0)
            u0 = 0;
        else
        {
            u0 = 2 - (l0 % 4);
            if ((l0 == 3 || l0 == 5) && l1 % 4 == 2)
                u0 = -u0;
        }
        if (l1 % 2 == 0)
            u1 = 0;
        else
        {
            u1 = 2 - (l1 % 4);
            if ((l1 == 3 || l1 == 5) && l0 % 4 == 2)
                u1 = -u1;
        }
#ifndef XMEN_ALWAYS_BINARY
        if (xmen_mip->base == xmen_mip->base2)
        {
#endif
            if (u0 > 0)
                xmen_addbit(u0p, j);
            if (u0 < 0)
                xmen_addbit(u0m, j);
            if (u1 > 0)
                xmen_addbit(u1p, j);
            if (u1 < 0)
                xmen_addbit(u1m, j);

#ifndef XMEN_ALWAYS_BINARY
        }
        else
        {
            if (u0 > 0)
                XMen_Add(u0p, xmen_mip->w1, u0p);
            if (u0 < 0)
                XMen_Add(u0m, xmen_mip->w1, u0m);
            if (u1 > 0)
                XMen_Add(u1p, xmen_mip->w1, u1p);
            if (u1 < 0)
                XMen_Add(u1m, xmen_mip->w1, u1m);
        }
#endif

        if (d0 + d0 == 1 + u0)
            d0 = 1 - d0;
        if (d1 + d1 == 1 + u1)
            d1 = 1 - d1;

        subdiv(xmen_mip->w2, 2, xmen_mip->w2);
        subdiv(xmen_mip->w3, 2, xmen_mip->w3);

#ifndef XMEN_ALWAYS_BINARY
        if (xmen_mip->base == xmen_mip->base2)
#endif
            j++;
#ifndef XMEN_ALWAYS_BINARY
        else
            premult(xmen_mip->w1, 2, xmen_mip->w1);
#endif
    }
    XMEN_OUT
    return;
}

/**
 * @description Performs elliptic curve point addition and subtraction between two points.
 *        This function computes both the sum and difference of two elliptic curve points `P` and `Q`,
 *        storing the results in `PP` (for addition) and `PM` (for subtraction). It handles special cases,
 *        such as when the points are equal or negatives of each other, and ensures that the resulting points
 *        are normalized. The function also validates the input points to ensure they are not in a general state.
 * @param[in] P   : Pointer to the first elliptic curve point for addition and subtraction.
 * @param[in] Q   : Pointer to the second elliptic curve point for addition and subtraction.
 * @param[out] PP : Pointer to the elliptic curve point where the result of `P + Q` will be stored.
 * @param[out] PM : Pointer to the elliptic curve point where the result of `P - Q` will be stored.
 * @return BOOL   : Returns TRUE if the operation was successful, otherwise returns FALSE.
 * @retval TRUE   : Both addition and subtraction were successfully computed.
 * @retval FALSE  : An error occurred, such as invalid parameters or encountering an exceptional condition.
 */
static BOOL ecurve_add_sub(XMen_ECurvePointT *P, XMen_ECurvePointT *Q, XMen_ECurvePointT *PP, XMen_ECurvePointT *PM)
{
    XMen_BigTPtr t1, t2, lam;

    if (xmen_mip->ERNUM)
        return FALSE;

    if (P->marker == XMEN_EPOINT_GENERAL || Q->marker == XMEN_EPOINT_GENERAL)
    {
        EasyXMen_berror(XMEN_ERR_BAD_PARAMETERS);
        XMEN_OUT
        return FALSE;
    }

    if (XMen_Compare(P->X, Q->X) == 0)
    { /* P=Q or P=-Q - shouldn't happen */
        epoint_copy(P, PP);
        ecurve_add(Q, PP);
        epoint_copy(P, PM);
        ecurve_sub(Q, PM);

        XMEN_OUT
        return TRUE;
    }

    t1 = xmen_mip->w10;
    t2 = xmen_mip->w11;
    lam = xmen_mip->w13;

    copy(P->X, t2);
    nres_modsub(t2, Q->X, t2);

    redc(t2, t2);
    invmodp(t2, xmen_mip->modulus, t2);
    nres(t2, t2);

    nres_modadd(P->X, Q->X, PP->X);
    copy(PP->X, PM->X);

    copy(P->Y, t1);
    nres_modsub(t1, Q->Y, t1);
    copy(t1, lam);
    nres_modmult(lam, t2, lam);
    copy(lam, t1);
    nres_modmult(t1, t1, t1);
    nres_modsub(t1, PP->X, PP->X);
    copy(Q->X, PP->Y);
    nres_modsub(PP->Y, PP->X, PP->Y);
    nres_modmult(PP->Y, lam, PP->Y);
    nres_modsub(PP->Y, Q->Y, PP->Y);

    copy(P->Y, t1);
    nres_modadd(t1, Q->Y, t1);
    copy(t1, lam);
    nres_modmult(lam, t2, lam);
    copy(lam, t1);
    nres_modmult(t1, t1, t1);
    nres_modsub(t1, PM->X, PM->X);
    copy(Q->X, PM->Y);
    nres_modsub(PM->Y, PM->X, PM->Y);
    nres_modmult(PM->Y, lam, PM->Y);
    nres_modadd(PM->Y, Q->Y, PM->Y);

    PP->marker = XMEN_EPOINT_NORMALIZED;
    PM->marker = XMEN_EPOINT_NORMALIZED;

    return TRUE;
}

/**
 * @description Handles big number error conditions.
 *        This function sets an error number in the global context (`xmen_mip`) when an error occurs,
 *        provided that error checking is enabled. It does not perform any action if error checking is disabled.
 * @param[in] nerr : The error code to set.
 * @return void
 */
static void EasyXMen_berror(int nerr)
{ /*  Big number error routine  */
    if (xmen_mip->ERCON)
    {
        xmen_mip->ERNUM = nerr;
        return;
    }
}

/**
 * @description Allocates and initializes memory for the first SM2 type structure.
 *        This function allocates memory for a single instance of the `XMen_SM2Type` structure and initializes it to zero.
 * @return XMen_SM2Type* : Pointer to the allocated and initialized `XMen_SM2Type` structure, or NULL on failure.
 */
static XMen_SM2Type *EasyXMen_first_alloc()
{
    return (XMen_SM2Type *)calloc(1, sizeof(XMen_SM2Type));
}

/**
 * @description Allocates and initializes memory for multi-precision operations.
 *        This function allocates a block of memory of specified size and initializes it to zero.
 *        If allocation fails, it records an out-of-memory error using the error handling routine.
 * @param[in] num : The number of elements to allocate.
 * @param[in] size: The size of each element in bytes.
 * @return void*  : Pointer to the allocated memory block, or NULL on failure.
 */
static void *EasyXMen_alloc(int num, int size)
{
    char *p;

    p = (char *)calloc(num, size);
    if (p == NULL)
        EasyXMen_berror(XMEN_ERR_OUT_OF_MEMORY);
    return (void *)p;
}

/**
 * @description Frees previously allocated memory.
 *        This function releases a block of memory that was previously allocated by `EasyXMen_alloc`.
 *        It performs no action if the provided address is NULL.
 * @param[in] addr : Pointer to the memory block to free.
 * @return void
 */
static void EasyXMen_free(void *addr)
{
    if (addr == NULL)
        return;
    free(addr);
    return;
}

/**
 * @description Allocates memory for multi-precision integer operations with appropriate reservation.
 *        This function calculates the necessary size for multi-precision integers based on the system's configuration
 *        and then allocates the required memory block using `EasyXMen_alloc`.
 * @param[in] num : The number of multi-precision integer elements to allocate.
 * @return void*  : Pointer to the allocated memory block, or NULL on failure.
 */
static void *EasyXMen_memalloc(int num)
{
    return EasyXMen_alloc(xmen_big_reserve(num, xmen_mip->nib - 1), 1);
}

/**
 * @description Sets the base for multi-precision arithmetic operations.
 *        This function configures the base used in multi-precision arithmetic, determining how many digits can fit into each computer word.
 *        It optimizes the base setting to pack as many digits as possible into each word while ensuring that calculations remain efficient and accurate.
 * @param[in] nb : The desired base value to set.
 * @return XMen_SmallType : Returns a status code or computed inverse base depending on the configuration.
 * @retval 0 : Base setting was successful or not applicable under certain configurations.
 * @retval xmen_mip->inverse_base : Inverse of the final base value when floating-point rounding is enabled.
 */
static XMen_SmallType EasyXMen_setbase(XMen_SmallType nb)
{ /* set base. Pack as many digits as  *
   * possible into each computer word  */
    XMen_SmallType temp;

#ifndef XMEN_NOFULLWIDTH
    BOOL fits;
    int bits;

    fits = FALSE;
    bits = MIRACL;
    while (bits > 1)
    {
        bits /= 2;
        temp = ((XMen_SmallType)1 << bits);
        if (temp == nb)
        {
            fits = TRUE;
            break;
        }
        if (temp < nb || (bits % 2) != 0)
            break;
    }
    if (fits)
    {
        xmen_mip->apbase = nb;
        xmen_mip->pack = MIRACL / bits;
        xmen_mip->base = 0;
        return 0;
    }
#endif
    xmen_mip->apbase = nb;
    xmen_mip->pack = 1;
    xmen_mip->base = nb;
#ifdef XMEN_SIMPLE_BASE
    return 0;
#else
    if (xmen_mip->base == 0)
        return 0;
    temp = XMEN_DIV(MAXBASE, nb);
    while (temp >= nb)
    {
        temp = XMEN_DIV(temp, nb);
        xmen_mip->base *= nb;
        xmen_mip->pack++;
    }
#ifdef XMEN_FP_ROUNDING
    xmen_mip->inverse_base = xmen_invert(xmen_mip->base);
    return xmen_mip->inverse_base;
#else
    return 0;
#endif
#endif
}

/**
 * @description Initializes the basic multi-precision arithmetic system configuration.
 *        This function sets up the multi-precision arithmetic environment by configuring various parameters,
 *        including the base, error checking, and workspace allocation. It also initializes random number generation
 *        and other critical components for subsequent operations.
 * @param[in,out] xmen_mip : Pointer to the multi-precision integer context (`XMen_SM2Type`) to initialize.
 * @param[in] nd         : The desired degree or size of multi-precision numbers.
 * @param[in] nb         : The base value to set for multi-precision arithmetic.
 * @return XMen_SM2Type* : Returns the initialized multi-precision integer context pointer.
 */
static XMen_SM2Type *EasyXMen_mirsys_basic(XMen_SM2Type *xmen_mip, int nd, XMen_SmallType nb)
{
#ifndef XMEN_NO_RAND
    int i;
#endif

    XMen_SmallType b, nw;

    if (xmen_mip == NULL)
        return NULL;

        /* digest hardware configuration */

#ifdef XMEN_NO_STANDARD_IO
    xmen_mip->ERCON = TRUE;
#else
    xmen_mip->ERCON = FALSE;
#endif
#ifndef XMEN_STATIC
    xmen_mip->logN = 0;
    xmen_mip->degree = 0;
    xmen_mip->chin.NP = 0;
#endif

    xmen_mip->user = NULL;
    xmen_mip->same = FALSE;
    xmen_mip->first_one = FALSE;
    xmen_mip->debug = FALSE;
    xmen_mip->AA = 0;
    xmen_mip->coord = 0;

#ifdef XMEN_NOFULLWIDTH
    if (nb == 0)
    {
        EasyXMen_berror(XMEN_ERR_BAD_BASE);
        XMEN_OUT
        return xmen_mip;
    }
#endif

#ifdef XMen_DltType
#ifndef XMEN_NOFULLWIDTH
    if (sizeof(XMen_DltType) < 2 * sizeof(int))
    { /* double length type, isn't */
        EasyXMen_berror(XMEN_ERR_NOT_DOUBLE_LEN);
        XMEN_OUT
        return xmen_mip;
    }
#endif
#endif

    if (nb == 1 || nb > MAXBASE)
    {
        EasyXMen_berror(XMEN_ERR_BAD_BASE);
        XMEN_OUT
        return xmen_mip;
    }

#ifdef XMEN_FP_ROUNDING
    if (xmen_setbase(nb) == 0)
    { /* unable in fact to control FP rounding */
        EasyXMen_berror(XMEN_ERR_NO_ROUNDING);
        XMEN_OUT
        return xmen_mip;
    }
#else
    EasyXMen_setbase(nb);
#endif

    b = xmen_mip->base;

#ifdef XMEN_SIMPLE_BASE
    if (b != 0)
    {
        EasyXMen_berror(XMEN_ERR_BAD_BASE);
        XMEN_OUT
        return xmen_mip;
    }
#endif

    xmen_mip->lg2b = 0;
    xmen_mip->base2 = 1;
#ifndef XMEN_SIMPLE_BASE
    if (b == 0)
    {
#endif
        xmen_mip->lg2b = MIRACL;
        xmen_mip->base2 = 0;
#ifndef XMEN_SIMPLE_BASE
    }
    else
        while (b > 1)
        {
            b = XMEN_DIV(b, 2);
            xmen_mip->lg2b++;
            xmen_mip->base2 *= 2;
        }
#endif

    xmen_mip->base_mask = xmen_mip->base - 1;

#ifdef XMEN_ALWAYS_BINARY
    if (xmen_mip->base != xmen_mip->base2)
    {
        EasyXMen_berror(XMEN_ERR_NOT_BINARY);
        XMEN_OUT
        return xmen_mip;
    }
#endif

    /* calculate total space for bigs */
    /*
     big -> |int len|small *ptr| alignment space | size in words +1| alignment up to multiple of 4 |
    */
    if (nd > 0)
        nw = XMEN_ROUNDUP(nd, xmen_mip->pack);
    else
        nw = XMEN_ROUNDUP(8 * (-nd), xmen_mip->lg2b);

    if (nw < 1)
        nw = 1;
    xmen_mip->nib = (int)(nw + 1); /* add one extra word for small overflows */

#ifdef XMEN_STATIC
    if (nw > XMEN_STATIC)
    {
        EasyXMen_berror(XMEN_ERR_TOO_BIG);
        XMEN_OUT
        return xmen_mip;
    }
#endif

    xmen_mip->nib = (int)(nw + 1); /* add one extra word for small overflows */

#ifdef XMEN_FLASH
    xmen_mip->workprec = xmen_mip->nib;
    xmen_mip->stprec = xmen_mip->nib;
    while (xmen_mip->stprec > 2 && xmen_mip->stprec > XMEN_FLASH / xmen_mip->lg2b)
        xmen_mip->stprec = (xmen_mip->stprec + 1) / 2;
    if (xmen_mip->stprec < 2)
        xmen_mip->stprec = 2;

#endif

#ifndef XMEN_DOUBLE_BIG
    xmen_mip->check = ON;
#else
    xmen_mip->check = OFF;
#endif

#ifndef XMEN_SIMPLE_BASE
#ifndef XMEN_SIMPLE_IO
    xmen_mip->IOBASE = 10; /* defaults */
#endif
#endif
    xmen_mip->ERNUM = 0;

    xmen_mip->NTRY = 6;
    xmen_mip->MONTY = ON;

#ifndef XMEN_SIMPLE_IO
    xmen_mip->INPLEN = 0;
    xmen_mip->IOBSIZ = XMEN_DEFAULT_BUFFER_SIZE;
#endif

#ifdef XMEN_STATIC
    xmen_mip->PRIMES = xmen_small_primes;
#else
    xmen_mip->PRIMES = NULL;
#ifndef XMEN_SIMPLE_IO
    xmen_mip->IOBUFF = (char *)EasyXMen_alloc(XMEN_DEFAULT_BUFFER_SIZE + 1, 1);
#endif
#endif
#ifndef XMEN_SIMPLE_IO
    xmen_mip->IOBUFF[0] = '\0';
#endif
    xmen_mip->qnr = 0;
    xmen_mip->cnr = 0;
    xmen_mip->TWIST = 0;
    xmen_mip->pmod8 = 0;
    xmen_mip->pmod9 = 0;

    /* quick start for rng. XMen_Irand(.) should be called first before serious use.. */

#ifndef XMEN_NO_RAND
    xmen_mip->ira[0] = 0x55555555;
    xmen_mip->ira[1] = 0x12345678;

    for (i = 2; i < NK; i++)
        xmen_mip->ira[i] = xmen_mip->ira[i - 1] + xmen_mip->ira[i - 2] + 0x1379BDF1;
    xmen_mip->rndptr = NK;
    xmen_mip->borrow = 0;
#endif

    xmen_mip->nib = 2 * xmen_mip->nib + 1;
#ifdef XMEN_FLASH
    if (xmen_mip->nib != (xmen_mip->nib & (XMEN_MSK)))
#else
    if (xmen_mip->nib != (int)(xmen_mip->nib & (XMEN_OBITS)))
#endif
    {
        EasyXMen_berror(XMEN_ERR_TOO_BIG);
        xmen_mip->nib = (xmen_mip->nib - 1) / 2;
        XMEN_OUT
        return xmen_mip;
    }
#ifndef XMEN_STATIC
    xmen_mip->workspace = (char *)EasyXMen_memalloc(XMEN_SPACES); /* grab workspace */
#else
    memset(xmen_mip->workspace, 0, XMEN_BIG_RESERVE(XMEN_SPACES));
#endif

    xmen_mip->M = 0;
    xmen_mip->fin = FALSE;
    xmen_mip->fout = FALSE;
    xmen_mip->active = ON;

    xmen_mip->nib = (xmen_mip->nib - 1) / 2;

    /* allocate memory for workspace variables */
#ifndef XMEN_DOUBLE_BIG

    xmen_mip->w0 = mirvar_mem(xmen_mip->workspace, 0); /* double length */
    xmen_mip->w1 = mirvar_mem(xmen_mip->workspace, 2);
    xmen_mip->w2 = mirvar_mem(xmen_mip->workspace, 3);
    xmen_mip->w3 = mirvar_mem(xmen_mip->workspace, 4);
    xmen_mip->w4 = mirvar_mem(xmen_mip->workspace, 5);
    xmen_mip->w5 = mirvar_mem(xmen_mip->workspace, 6);  /* double length */
    xmen_mip->w6 = mirvar_mem(xmen_mip->workspace, 8);  /* double length */
    xmen_mip->w7 = mirvar_mem(xmen_mip->workspace, 10); /* double length */
    xmen_mip->w8 = mirvar_mem(xmen_mip->workspace, 12);
    xmen_mip->w9 = mirvar_mem(xmen_mip->workspace, 13);
    xmen_mip->w10 = mirvar_mem(xmen_mip->workspace, 14);
    xmen_mip->w11 = mirvar_mem(xmen_mip->workspace, 15);
    xmen_mip->w12 = mirvar_mem(xmen_mip->workspace, 16);
    xmen_mip->w13 = mirvar_mem(xmen_mip->workspace, 17);
    xmen_mip->w14 = mirvar_mem(xmen_mip->workspace, 18);
    xmen_mip->w15 = mirvar_mem(xmen_mip->workspace, 19);
    xmen_mip->sru = mirvar_mem(xmen_mip->workspace, 20);
    xmen_mip->modulus = mirvar_mem(xmen_mip->workspace, 21);
    xmen_mip->pR = mirvar_mem(xmen_mip->workspace, 22); /* double length */
    xmen_mip->A = mirvar_mem(xmen_mip->workspace, 24);
    xmen_mip->B = mirvar_mem(xmen_mip->workspace, 25);
    xmen_mip->one = mirvar_mem(xmen_mip->workspace, 26);
#ifdef XMEN_FLASH
    xmen_mip->pi = mirvar_mem(xmen_mip->workspace, 27);
#endif

#else
    /* w0-w7 are double normal length */
    xmen_mip->w0 = mirvar_mem(xmen_mip->workspace, 0); /* quad length */
    xmen_mip->w1 = mirvar_mem(xmen_mip->workspace, 4); /* double length */
    xmen_mip->w2 = mirvar_mem(xmen_mip->workspace, 6);
    xmen_mip->w3 = mirvar_mem(xmen_mip->workspace, 8);
    xmen_mip->w4 = mirvar_mem(xmen_mip->workspace, 10);
    xmen_mip->w5 = mirvar_mem(xmen_mip->workspace, 12); /* quad length */
    xmen_mip->w6 = mirvar_mem(xmen_mip->workspace, 16); /* quad length */
    xmen_mip->w7 = mirvar_mem(xmen_mip->workspace, 20); /* quad length */
    xmen_mip->w8 = mirvar_mem(xmen_mip->workspace, 24);

    xmen_mip->w9 = mirvar_mem(xmen_mip->workspace, 25);
    xmen_mip->w10 = mirvar_mem(xmen_mip->workspace, 26);
    xmen_mip->w11 = mirvar_mem(xmen_mip->workspace, 27);
    xmen_mip->w12 = mirvar_mem(xmen_mip->workspace, 28);
    xmen_mip->w13 = mirvar_mem(xmen_mip->workspace, 29);
    xmen_mip->w14 = mirvar_mem(xmen_mip->workspace, 30);
    xmen_mip->w15 = mirvar_mem(xmen_mip->workspace, 31);
    xmen_mip->sru = mirvar_mem(xmen_mip->workspace, 32);
    xmen_mip->modulus = mirvar_mem(xmen_mip->workspace, 33);
    xmen_mip->pR = mirvar_mem(xmen_mip->workspace, 34); /* double length */
    xmen_mip->A = mirvar_mem(xmen_mip->workspace, 36);
    xmen_mip->B = mirvar_mem(xmen_mip->workspace, 37);
    xmen_mip->one = mirvar_mem(xmen_mip->workspace, 38);
#ifdef XMEN_FLASH
    xmen_mip->pi = mirvar_mem(xmen_mip->workspace, 39);
#endif

#endif
    XMEN_OUT
    return xmen_mip;
}

/**
 * @description Sets a big or flash number to zero.
 *        This function initializes all elements of the given multi-precision number (`flash`) to zero and sets its length to zero.
 * @param[in,out] x : Pointer to the multi-precision number to be set to zero.
 * @return void
 */
static void zero(flash x)
{ /* set big/flash number to zero */
    int i, n;
    XMen_SmallType *g;
    if (x == NULL)
        return;
#ifdef XMEN_FLASH
    n = xmen_lent(x);
#else
    n = (x->len & XMEN_OBITS);
#endif
    g = x->w;

    for (i = 0; i < n; i++)
        g[i] = 0;

    x->len = 0;
}

/**
 * @description Checks if the provided multi-precision number is a Flash type.
 *        This function determines whether the given multi-precision number (`flash`) is of the Flash type, which allows for more efficient arithmetic operations.
 * @param[in] x : Pointer to the multi-precision number to check.
 * @return BOOL : Returns TRUE if `x` is a Flash type; otherwise, returns FALSE.
 */
static BOOL xmen_notint(flash x)
{ /* returns TRUE if x is Flash */
#ifdef XMEN_FLASH
    if ((((x->len & (XMEN_OBITS)) >> (XMEN_BTS)) & (XMEN_MSK)) != 0)
        return TRUE;
#endif
    return FALSE;
}

/**
 * @description Copies the contents of one multi-precision number to another.
 *        This function copies the value of the source multi-precision number (`x`) to the destination multi-precision number (`y`).
 *        It ensures that the destination is properly initialized and handles cases where either input might be NULL.
 * @param[in] x : Pointer to the source multi-precision number.
 * @param[in,out] y : Pointer to the destination multi-precision number.
 * @return void
 */
static void copy(flash x, flash y)
{ /* copy x to y: y=x  */
    int i, nx, ny;
    XMen_SmallType *gx, *gy;
    if (x == y || y == NULL)
        return;

    if (x == NULL)
    {
        zero(y);
        return;
    }

#ifdef XMEN_FLASH
    ny = xmen_lent(y);
    nx = xmen_lent(x);
#else
    ny = (y->len & (XMEN_OBITS));
    nx = (x->len & (XMEN_OBITS));
#endif

    gx = x->w;
    gy = y->w;

    for (i = nx; i < ny; i++)
        gy[i] = 0;
    for (i = 0; i < nx; i++)
        gy[i] = gx[i];
    y->len = x->len;
}

/**
 * @description Performs multiplication and division on small types.
 *        This function multiplies two small type values (`a` and `b`), adds a third value (`c`), and returns the high part of the result,
 *        storing the low part in the provided pointer (`rp`).
 * @param[in] a : First multiplier.
 * @param[in] b : Second multiplier.
 * @param[in] c : Value to add after multiplication.
 * @param[out] rp : Pointer to store the low part of the result.
 * @return XMen_SmallType : The high part of the result from the multiplication and addition.
 */
static XMen_SmallType muldvd(XMen_SmallType a, XMen_SmallType b, XMen_SmallType c, XMen_SmallType *rp)
{
    union doubleword dble;
    dble.d = (XMen_LargeType)a * b + c;

    *rp = dble.h[BYTE_BOT];
    return dble.h[BYTE_TOP];
}

/**
 * @description Performs multiplication, division, and modulus operations on small types.
 *        This function multiplies two small type values (`a` and `b`), adds a third value (`c`), divides the result by `m`,
 *        and returns the quotient while storing the remainder in the provided pointer (`rp`).
 * @param[in] a : First multiplier.
 * @param[in] b : Second multiplier.
 * @param[in] c : Value to add after multiplication.
 * @param[in] m : Divisor for the division operation.
 * @param[out] rp : Pointer to store the remainder of the division.
 * @return XMen_SmallType : The quotient from the division operation.
 */
static XMen_SmallType muldiv(XMen_SmallType a, XMen_SmallType b, XMen_SmallType c, XMen_SmallType m, XMen_SmallType *rp)
{
    XMen_SmallType q;
    // XMen_LargeType ldres;
    p = (XMen_LargeType)a * b + c;
    q = (XMen_SmallType)(XMEN_LROUND(p / m));
    *rp = (XMen_SmallType)(p - (XMen_LargeType)q * m);
    return q;
}

/**
 * @description Multiplies a multi-precision number by a small integer.
 *        This function multiplies the multi-precision number `x` by a small integer `sn` and stores the result in `z`.
 *        It handles cases where `x` and `z` are the same, ensuring proper initialization and overflow checking.
 * @param[in] x : Pointer to the source multi-precision number to be multiplied.
 * @param[in] sn : The small integer multiplier.
 * @param[out] z : Pointer to the destination multi-precision number where the result will be stored.
 * @return void
 */
static void xmen_pmul(XMen_BigTPtr x, XMen_SmallType sn, XMen_BigTPtr z)
{
    int m, xl;
    unsigned int sx;
    XMen_SmallType carry, *xg, *zg;

#ifdef XMEN_ITANIUM
    XMen_SmallType tm;
#endif
#ifdef XMEN_WIN64
    XMen_SmallType tm;
#endif
    union doubleword dble;
    XMen_LargeType dbled;
    XMen_LargeType ldres;

    if (x != z)
    {
        zero(z);
        if (sn == 0)
            return;
    }
    else if (sn == 0)
    {
        zero(z);
        return;
    }
    m = 0;
    carry = 0;
    sx = x->len & XMEN_MSBIT;
    xl = (int)(x->len & XMEN_OBITS);

#ifndef XMEN_SIMPLE_BASE
    if (xmen_mip->base == 0)
    {
#endif
#ifndef XMEN_NOFULLWIDTH
        xg = x->w;
        zg = z->w;
        for (m = 0; m < xl; m++)
        {
            dble.d = (XMen_LargeType)x->w[m] * sn + carry;
            carry = dble.h[BYTE_TOP];
            z->w[m] = dble.h[BYTE_BOT];
        }
        if (carry > 0)
        {
            m = xl;
            if (m >= xmen_mip->nib && xmen_mip->check)
            {
                EasyXMen_berror(XMEN_ERR_OVERFLOW);
                return;
            }
            z->w[m] = carry;
            z->len = m + 1;
        }
        else
            z->len = xl;
#endif
#ifndef XMEN_SIMPLE_BASE
    }
    else
        while (m < xl || carry > 0)
        { /* multiply each digit of x by n */

            if (m > xmen_mip->nib && xmen_mip->check)
            {
                EasyXMen_berror(XMEN_ERR_OVERFLOW);
                return;
            }
            dbled = (XMen_LargeType)x->w[m] * sn + carry;
#ifdef XMEN_FP_ROUNDING
            carry = (XMen_SmallType)XMEN_LROUND(dbled * xmen_mip->inverse_base);
#else
            if (xmen_mip->base == xmen_mip->base2)
                carry = (XMen_SmallType)(dbled >> xmen_mip->lg2b);
            else
                carry = (XMen_SmallType)XMEN_LROUND(dbled / xmen_mip->base);
#endif
            z->w[m] = (XMen_SmallType)(dbled - (XMen_LargeType)carry * xmen_mip->base);
            m++;
            z->len = m;
        }
#endif
    if (z->len != 0)
        z->len |= sx;
}

/**
 * @description Converts an unsigned integer to a multi-precision number format.
 *        This function takes an unsigned integer `n` and converts it into the multi-precision number format (`XMen_BigTPtr`),
 *        storing the result in `x`. It handles different base configurations and ensures proper conversion for various cases.
 * @param[in] n : The unsigned integer to convert.
 * @param[out] x : Pointer to the multi-precision number where the result will be stored.
 * @return void
 */
static void uconvert(unsigned int n, XMen_BigTPtr x)
{ /*  convert unsigned integer n to big number format  */
    int m;

    zero(x);
    if (n == 0)
        return;

    m = 0;
#ifndef XMEN_SIMPLE_BASE
    if (xmen_mip->base == 0)
    {
#endif
#ifndef XMEN_NOFULLWIDTH
#if XMEN_IBITS > MIRACL
        while (n > 0)
        {
            x->w[m++] = (XMen_SmallType)(n % ((XMen_SmallType)1 << (MIRACL)));
            n /= ((XMen_SmallType)1 << (MIRACL));
        }
#else
        x->w[m++] = (XMen_SmallType)n;
#endif
#endif
#ifndef XMEN_SIMPLE_BASE
    }
    else
        while (n > 0)
        {
            x->w[m++] = XMEN_REMAIN((XMen_SmallType)n, xmen_mip->base);
            n /= xmen_mip->base;
        }
#endif
    x->len = m;
}

static int exsign(flash x)
{ /* extract sign of big/flash number */
    if ((x->len & (XMEN_MSBIT)) == 0)
        return PLUS;
    else
        return MINUS;
}

/**
 * @description Adds two multi-precision numbers.
 *        This function adds two positive multi-precision numbers (`x` and `y`) and stores the result in `z`.
 *        It handles cases where `x`, `y`, and `z` might point to the same memory location and ensures proper overflow checking.
 * @param[in] x : Pointer to the first multi-precision number to add.
 * @param[in] y : Pointer to the second multi-precision number to add.
 * @param[out] z : Pointer to the multi-precision number where the sum will be stored.
 * @return void
 */
static void xmen_padd(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr z)
{ /*  add two  big numbers, z=x+y where *
   *  x and y are positive              */
    int i, lx, ly, lz, la;
    XMen_SmallType carry, psum;
    XMen_SmallType *gx, *gy, *gz;

    lx = (int)x->len;
    ly = (int)y->len;

    if (ly > lx)
    {
        lz = ly;
        la = lx;
        if (x != z)
            copy(y, z);
        else
            la = ly;
    }
    else
    {
        lz = lx;
        la = ly;
        if (y != z)
            copy(x, z);
        else
            la = lx;
    }
    carry = 0;
    z->len = lz;
    gx = x->w;
    gy = y->w;
    gz = z->w;
    if (lz < xmen_mip->nib || !xmen_mip->check)
        z->len++;
#ifndef XMEN_SIMPLE_BASE
    if (xmen_mip->base == 0)
    {
#endif
        for (i = 0; i < la; i++)
        { /* add by columns to length of the smaller number */
            psum = gx[i] + gy[i] + carry;
            if (psum > gx[i])
                carry = 0;
            else if (psum < gx[i])
                carry = 1;
            gz[i] = psum;
        }
        for (; i < lz && carry > 0; i++)
        { /* add by columns to the length of larger number (if there is a carry) */
            psum = gx[i] + gy[i] + carry;
            if (psum > gx[i])
                carry = 0;
            else if (psum < gx[i])
                carry = 1;
            gz[i] = psum;
        }
        if (carry)
        { /* carry left over - possible overflow */
            if (xmen_mip->check && i >= xmen_mip->nib)
            {
                EasyXMen_berror(XMEN_ERR_OVERFLOW);
                return;
            }
            gz[i] = carry;
        }
#ifndef XMEN_SIMPLE_BASE
    }
    else
    {
        for (i = 0; i < la; i++)
        { /* add by columns */
            psum = gx[i] + gy[i] + carry;
            carry = 0;
            if (psum >= xmen_mip->base)
            { /* set carry */
                carry = 1;
                psum -= xmen_mip->base;
            }
            gz[i] = psum;
        }
        for (; i < lz && carry > 0; i++)
        {
            psum = gx[i] + gy[i] + carry;
            carry = 0;
            if (psum >= xmen_mip->base)
            { /* set carry */
                carry = 1;
                psum -= xmen_mip->base;
            }
            gz[i] = psum;
        }
        if (carry)
        { /* carry left over - possible overflow */
            if (xmen_mip->check && i >= xmen_mip->nib)
            {
                EasyXMen_berror(XMEN_ERR_OVERFLOW);
                return;
            }
            gz[i] = carry;
        }
    }
#endif
    if (gz[z->len - 1] == 0)
        z->len--;
}

/**
 * @description Subtracts two multi-precision numbers.
 *        This function subtracts the second multi-precision number (`y`) from the first one (`x`) and stores the result in `z`.
 *        It assumes that both `x` and `y` are positive, and `x` is greater than `y`. If these conditions are not met, it reports an error.
 * @param[in] x : Pointer to the first multi-precision number (minuend).
 * @param[in] y : Pointer to the second multi-precision number (subtrahend).
 * @param[out] z : Pointer to the multi-precision number where the difference will be stored.
 * @return void
 */
static void xmen_psub(XMen_BigTPtr x, XMen_BigTPtr y, XMen_BigTPtr z)
{ /*  subtract two big numbers z=x-y      *
   *  where x and y are positive and x>y  */
    int i, lx, ly;
    XMen_SmallType borrow, pdiff;
    XMen_SmallType *gx, *gy, *gz;

    lx = (int)x->len;
    ly = (int)y->len;
    if (ly > lx)
    {
        EasyXMen_berror(XMEN_ERR_NEG_RESULT);
        return;
    }
    if (y != z)
        copy(x, z);
    else
        ly = lx;
    z->len = lx;
    gx = x->w;
    gy = y->w;
    gz = z->w;
    borrow = 0;
#ifndef XMEN_SIMPLE_BASE
    if (xmen_mip->base == 0)
    {
#endif
        for (i = 0; i < ly || borrow > 0; i++)
        { /* subtract by columns */
            if (i > lx)
            {
                EasyXMen_berror(XMEN_ERR_NEG_RESULT);
                return;
            }
            pdiff = gx[i] - gy[i] - borrow;
            if (pdiff < gx[i])
                borrow = 0;
            else if (pdiff > gx[i])
                borrow = 1;
            gz[i] = pdiff;
        }
#ifndef XMEN_SIMPLE_BASE
    }
    else
        for (i = 0; i < ly || borrow > 0; i++)
        { /* subtract by columns */
            if (i > lx)
            {
                EasyXMen_berror(XMEN_ERR_NEG_RESULT);
                return;
            }
            pdiff = gy[i] + borrow;
            borrow = 0;
            if (gx[i] >= pdiff)
                pdiff = gx[i] - pdiff;
            else
            { /* set borrow */
                pdiff = xmen_mip->base + gx[i] - pdiff;
                borrow = 1;
            }
            gz[i] = pdiff;
        }
#endif
    xmen_lzero(z);
}

/**
 * @description Performs the required addition or subtraction operation based on the signs of the operands.
 *        This function selects and executes either an addition or subtraction operation between two multi-precision numbers (`x` and `y`),
 *        storing the result in `z`. The operation performed depends on the values of `d` and the signs of `x` and `y`.
 * @param[in] x : Pointer to the first multi-precision number operand.
 * @param[in] d : Decision parameter that influences whether addition or subtraction is performed.
 * @param[in] y : Pointer to the second multi-precision number operand.
 * @param[out] z : Pointer to the multi-precision number where the result will be stored.
 * @return void
 */
static void xmen_select(XMen_BigTPtr x, int d, XMen_BigTPtr y, XMen_BigTPtr z)
{ /* perform required add or subtract operation */
    int sx, sy, sz, jf, xgty;
#ifdef XMEN_FLASH
    if (xmen_notint(x) || xmen_notint(y))
    {
        EasyXMen_berror(XMEN_ERR_INT_OP);
        return;
    }
#endif
    sx = exsign(x);
    sy = exsign(y);
    sz = 0;
    x->len &= XMEN_OBITS; /* force operands to be positive */
    y->len &= XMEN_OBITS;
    xgty = XMen_Compare(x, y);
    jf = (1 + sx) + (1 + d * sy) / 2;
    switch (jf)
    { /* branch according to signs of operands */
    case 0:
        if (xgty >= 0)
            xmen_padd(x, y, z);
        else
            xmen_padd(y, x, z);
        sz = MINUS;
        break;
    case 1:
        if (xgty <= 0)
        {
            xmen_psub(y, x, z);
            sz = PLUS;
        }
        else
        {
            xmen_psub(x, y, z);
            sz = MINUS;
        }
        break;
    case 2:
        if (xgty >= 0)
        {
            xmen_psub(x, y, z);
            sz = PLUS;
        }
        else
        {
            xmen_psub(y, x, z);
            sz = MINUS;
        }
        break;
    case 3:
        if (xgty >= 0)
            xmen_padd(x, y, z);
        else
            xmen_padd(y, x, z);
        sz = PLUS;
        break;
    }
    if (sz < 0)
        z->len ^= XMEN_MSBIT; /* set sign of result         */
    if (x != z && sx < 0)
        x->len ^= XMEN_MSBIT; /* restore signs to operands  */
    if (y != z && y != x && sy < 0)
        y->len ^= XMEN_MSBIT;
}

/**
 * @description Premultiplies a multi-precision number by an integer.
 *        This function multiplies the multi-precision number `x` by an integer `n` and stores the result in `z`.
 *        It handles special cases such as `n` being zero or one, and manages the sign of the result when `n` is negative.
 * @param[in] x : Pointer to the multi-precision number to be multiplied.
 * @param[in] n : The integer multiplier.
 * @param[out] z : Pointer to the multi-precision number where the result will be stored.
 * @return void
 */
static void premult(XMen_BigTPtr x, int n, XMen_BigTPtr z)
{ /* premultiply a big number by an int z=x.n */

    if (xmen_mip->ERNUM)
        return;

    XMEN_IN(9)

#ifdef XMEN_FLASH
    if (xmen_notint(x))
    {
        EasyXMen_berror(XMEN_ERR_INT_OP);
        XMEN_OUT
        return;
    }
#endif
    if (n == 0) /* test for some special cases  */
    {
        zero(z);
        XMEN_OUT
        return;
    }
    if (n == 1)
    {
        copy(x, z);
        XMEN_OUT
        return;
    }
    if (n < 0)
    {
        n = (-n);
        xmen_pmul(x, (XMen_SmallType)n, z);
        if (z->len != 0)
            z->len ^= XMEN_MSBIT;
    }
    else
        xmen_pmul(x, (XMen_SmallType)n, z);
    XMEN_OUT
}

/**
 * @description Extracts the numerator of a multi-precision number.
 *        This function extracts the numerator from the multi-precision number `x` and stores it in `y`.
 *        It handles cases where `x` is not an integer type and ensures proper handling of signs and special cases.
 * @param[in] x : Pointer to the multi-precision number from which the numerator will be extracted.
 * @param[out] y : Pointer to the multi-precision number where the numerator will be stored.
 * @return void
 */
static void numer(flash x, XMen_BigTPtr y)
{ /* extract numerator of x */
    int i, ln, ld;
    unsigned int s, ly;

    if (xmen_mip->ERNUM)
        return;
    if (xmen_notint(x))
    {
        s = (x->len & XMEN_MSBIT);
        ly = (x->len & XMEN_OBITS);
        ln = (int)(ly & XMEN_MSK);
        if (ln == 0)
        {
            if (s == XMEN_MSBIT)
                convert((-1), y);
            else
                convert(1, y);
            return;
        }
        ld = (int)((ly >> XMEN_BTS) & XMEN_MSK);
        if (x != y)
        {
            for (i = 0; i < ln; i++)
                y->w[i] = x->w[i];
            for (i = ln; i < xmen_lent(y); i++)
                y->w[i] = 0;
        }
        else
            for (i = 0; i < ld; i++)
                y->w[ln + i] = 0;
        y->len = (ln | s);
    }
    else
        copy(x, y);
}
